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A doubly curved element for a shell of revolution which has arbitrar- 
ily curved meridians is developed and analyzed. Meridional curvature is 
calculated using a highly accurate polynomial approximation. The dis- 
placement functions selected satisfy interelement compatibility and 
contain all the lower modes of a complete set of straining functions. 
Non-straining modes corresponding to rigid body motions are introduced 
into the final stiffness matrix for any conceivable rigid body motions. 

The direct stiffness method was used to construct a stress and 
strain analysis program and the results of the analysis of a few problems 
are compared to classical solutions to establish the integrity of the 



element. 
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I. INTRODUCTION 



A. GENERAL DESCRIPTION 

The concept of dividing complex structures into an assemblage of 
individual structural components or elements is axiomatic to the struc- 
tural engineer. However, the ability to execute a breaking down and 
reassembling process has been dependent upon the availability of an- 
alytical tools. Not until the present decade were developments in high- 
speed digital computers sufficient to make possible the fundamental 
approach to problems of structural analysis known as the finite element 
method . 

The finite element method is essentially the analysis of an idealized 
structural system composed of a finite number of elements interconnected 
at a finite number of joints or nodal points. It is the finite characteristic 
of structural connectivity which distinguishes the finite element problem 
from problems of continuum mechanics. 

B. HISTORICAL BACKGROUND 

The finite element method of analysis is applicable to one, two and 
three dimensional structures. Background discussion will be limited to 
the two dimensional case, although recently there has been considerable 
research activity on three dimensional structures . 

Various shapes of finite elements have been employed in plane stress 
and plane strain analysis. Early research was focused on the development 
of flat rectangular and triangular finite elements . An abundance of papers. 
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technical reports, and theses were produced; a complete bibliography 
would be beyond the scope of this paper. It is sufficient to say however, 
that requirements for efficient flat plate finite elements have been formu- 
lated and tested. Considerable success was reported in Refs. 1 and 2 
using these elements to approximate shells of free form. 

The development of curved finite elements , a considerably more 
difficult problem, has progressed at a much slower pace. Some early 
success was achieved in the development of truncated conical ring 
elements. However, problem formulation using these elements has not 
been generalized. Truncated conical elements have produced excellent 
results for particular problems but have also resulted in unsatisfactory 
analysis of similar problems. 

A ring element with arbitrarily curved meridians developed by R. E. 
Jones and D. R. Strome [Ref. 3] produced acceptable results for the case 
of shells of revolution under pressure loads. 

In 1968, G. Cantin [Ref. 4] successfully developed a curved finite 
element for cylindrical shells. This curved element was developed using 
6 nodal coordinates per joint (uj, Vj , Wj , w^^ , and implicitly 

contained a complete description of rigid body motions. This element was 
tested and found to produce excellent results for a broad class of cylin- 
drical shell and flat plate problems for both membrane and bending 
behavior. (The flat element is obtained by letting R-*”.) 
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C . OBJECTIVE 



The objective of the author's research is to develop the stiffness 
matrix of a four nodal point finite element for a shell of revolution which 
has arbitrarily curved meridians, hereafter to be called a curvilinear shell 
finite element. 
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II. THE FINITE ELEMENT METHOD 



A. GENERAL 

Finite element analysis of an elastic continuum logically separates 
into four phases: (1) structural idealization, (2) evaluation of the element 
stiffness properties, (3) the evaluation of the force -displacement be- 
havior of the overall element assemblage, and (4) the evaluation of 
element stress and strain. 

B. STRUCTURAL IDEALIZATION 

Structural idealization is the division of the actual elastic continuum 
into an assemblage of geometrically compatible structural elements which 
are interconnected at a discrete number of nodal points, through which 
the forces are transmitted. 

The discretization of the continuum is a physical process, subject 
only to conditions of geometric compatibility, rather than an approxima- 
tion of a mathematical nature. The idealized structure, then, is con- 
structed of elements which possess the material properties of the actual 
continuum . 

A basic reduction is necessary to make the idealized structure 
mathematically tractable. This reduction, based on an a priori selection 
of possible deformations of the elastic continuum, constrains the 
elements to deform in only a certain number of predictable shapes . These 
constraints called displacement functions uniquely define displacements 
and the state of strain within the element in terms of nodal displacements. 

These functions are the keystone of a successful idealization. 
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C . EVALUATION OF ELEMENT STIFFNESS PROPERTIES 

The evaluation of the stiffness properties of the element is the crit- 
ical phase of the finite element method. As the evaluation of the stiffness 
properties of a curvilinear shell finite element is the object of the 
author's research, procedures will be discussed in following sections. 

For the present, it is sufficient to say that the element stiffness matrix 
relates the nodal forces and nodal displacements . This relation takes 
the form 

{F]=[k]{6} (2.1) 

where 

(F] = nodal force vector 

{6} = nodal point displacement vector 

[k] = element stiffness matrix 

D. EVALUATION OF THE FORCE -DISPLACEMENT BEHAVIOR OF THE 
IDEALIZED STRUCTURE 

The force-displacement behavior of the overall element assemblage 
is obtained by the systematic superposition of individual element stiff- 
ness matrices and load vectors. This technique, called the "direct 
stiffness method," is well documented in the literature. 

E. EVALUATION OF ELEMENT STRESS AND STRAIN 

When nodal displacements are known^ element strain is directly 
obtained using classical shell theory strain-displacement relationships . 
Element stress is then calculated using constitutive laws. 
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III. CONSTRUCTION OF THE STIFFNESS MATRIX 



A. GENERAL 

The stiffness matrix for a particular finite elment is the relationship 
between a set of generalized displacements and the corresponding 
generalized forces. 

Standard matrix notation, square brackets denoting rectangular 
matrices, curly brackets representing column matrices, superscript (-1) 
denoting inversion, and superscript T denoting transposition, will be 
used . 

B. GEOMETRIC DERIVATION 

The geometric description of the curvilinear shell finite element 
must contain sufficient parameters to guarantee that an assemblage of 
elements will form a smooth surface whose slope is continuous and whose 
dimensions at element boundaries match those of the actual continuum. 
Furthermore, the geometric description must be sufficiently general to 
describe shells of revolution which have meridional inflection or reverse 
curvature as well as shells of degenerate geometry such as cylindrical 
shells and flat plates. An extension of the descriptive procedure used 
in Ref. 3 was found to satisfy the above requirements. 

Geometrical symbols for the curvilinear shell finite element are 
defined in Figure 1. Nondimensional meridional arc length is denoted 
by 
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NOTE- ALL DIMENSIONS BASED 
ON SHELL MIDSURFACE 

Fig. I DEFINITION OF GEOMETRICAL SYMBOLS 
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4= - 

^ z 



and nondimensional parallel arc length is denoted by 



(3.1) 



i 

oc 



(3.2) 



Through the use of a second order polynomial in 4 for the cosine of the 
angle 0 , 

cos 0 = Cq + Ci4 + (3.3) 

and the values of 

r(^), r^^^ , 0 , 0^^^, O' (see Figure 1) 

expressions for completely defining element geometry were directly 
formulated. The radius of the element is given by 



_ Ji) 



.4 



r(4) = fVi; + X j cos 0 d4 = r(0) + -^ (c^ 4 + i Ci4 + y C24®) 

(3.4) 

and the principal radii of curvature are 



d0/d4 



Z sin 0 

Cj + 



R2 = 



. . r .(4). . 

sin 0 



where 



c = cos 0 



Cl = 2 [3 



(i) 



,(2)_ ,(i) 



- 2 cos 0(^) - cos 0 ^ ] 



= 3 [-2 7^"’ )* ^ 



(3.5) 



(3.6) 



(3.7) 

(3.8) 

(3.9) 



14 



It is worthy of note that all values required for these expressions 
are found from engineering drawings or models. 

For particular cases, such as conical shells and cylindrical shells 
where geometry is defined by analytical expressions, geometric 
approximation is not required. 



C. DISPLACEMENT FUNCTIONS' 



The displacement functions 



{u} 



/ \ 
^ .V) ( 



\w(4 .V)J 



(3.10) 



must satisfy three critical requirements . 

1 . They must insure that continuity of the structure can be auto- 
matically maintained across common boundaries of adjacent elements. 

2 . They must insure that all the lower modes of a complete set of 
straining modes will be present in the stiffness matrix. 

3. They must contain an exact description of all possible rigid 
body motions of the element. 

It can be shown [Ref. 5] that satisfaction of the first two requirements 
will guarantee that the strain energy of the "assemblage of elements" 
approximation will represent, for specified loads,a lower bound to the 
strain energy of the actual continuum. Satisfaction of the third 
requirement is necessary to expose hidden constraints opposing rigid 
body motion. 
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1 . Application to the Curvilinear Shell Finite Element 



The polynomial expansion: 
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satisfies the requirements of interelement compatibility and completeness 
of strain field. Requirement three, that an exact description of all 
possible rigid body motions of the element be represented, is not 
satisfied, however, as will be shown in the following section. 

D. RIGID BODY MODES 

It has been established in Ref. 6, that conditions of equilibrium 
are not satisfied unless the displacement functions contain an exact 
description of rigid body motions of the element. Moreover, it has been 
established in Ref. 7 that inclusion of rigid body motion is essential and 
cannot be represented by independent displacement components. Rigid 
body motion, then, must be introduced without compromise to interelement 
compatibility requirements . 



16 



A general method for including rigid body motion is developed in 



Ref. 6. 

1 . Application to the Curvilinear Shell Finite Element 

Consider a point P, specified by position vector p, on the 
element in the system of reference shown in Figure 2 . If the element is 
submitted to general translation then, 

6^ = 6j + 63]' + 63k (3.13) 

in the system of reference (x, y, z) . The corresponding scalar compo- 
nents of the displacement field are: 
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- sin 0 
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Sin 0 cos 


6 


sin 0 
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cos 0 







or more compactly 

{u} = tl^TRANS^ 

When submitted to rigid body rotation of small amplitude 

/8 = /3ii + jSsJ + jSak (3.15) 

the displacement vector field of the element in reference system (x, y, z) 
is 

= iSx p (3.16) 



The corresponding scalar expression is 
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(3.17) 
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Fig 2 DEFINITION OF COORDINATE SYSTEMS 
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or more compactly 

{u} = [fill 

Combining expressions (3.14) and (3.17) gives: 



or 
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fu.) = [R] tUj" 1 

Matrix [R]is given explicitly in Table 1. 



(3.18) 



(3.19) 



E. CLASSICAL SHELL THEORY 

In the introduction it was emphasized that the difference between 
finite element problems and problems in continuum mechanics is based 
on a physical approximation only. Classical theories for shells from 
continuum mechanics are, therefore, applicable in the finite element 
approach . 

The consitutive law for an isotropic, homogenous, thin elastic shell 
can be written: 
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(3.20) 
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where 



Et 

;Ks = V Ki ; K3 = i Ki 



Ki = 

Dj^ — 22''’- ~ V T)i ; D3 = -5 (l-i/) D]^ 

or more compactly 

{N} = [E] {e} (3.21) 

Young's modulus E and Poisson's ratio u, according to definition 
are constants. The thickness t, however, is allowed to vary in the 
meridional direction. In this thesis, shells of uniform thickness and 
shells of linearly varying thickness are considered.^ 

Numerous theories exist for curvilinear shells. The venerable 
theories have been subject to years of critical inspection; their merits 
and limitations are well known. Several theories investigated would 
have been acceptable. Kraus' formulation [Ref. 8] was selected for the 
following reasons: 

1. The theory is consistent. That is, the theory holds for the 
limiting cases of cylindrical shell and flat plate elements. This quality 
is of basic importance since mesh size reduction to the extent where the 
element is nearly flat is necessary in many finite element problems. 

2 . The theory is presented in a form that can be directly applied 
to this study. 

3 . The theory allows strain free modes for rigid body motions . 

^ Any thickness which varies as a function of ^ could be considered 
by changing only one card in the computer program. 
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The strain displacement relations are: 
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(3.22) 



or more compactly 

{c}=[D]Cu} (3.23) 

All classical theories investigated by the author contained inherent 
singularities. Kraus' formulation was no exception. It can be seen that 
several terms in equation (3.22) become infinite whenever r is equal to 
zero. This theoretical inadequacy has been avoided by specifying either 
a small hole or a small rigid insert with appropriate boundary conditions , 
at the pole point. The problem solved in Chapter 5, Section 2 
demonstrates this method. 
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F . NODAL COORDINATES 



For computational convenience, physical nodal coordinates are 
transformed to a generalized coordinate system. This transformation 
facilitates the expression of displacements, elastic characteristics, and 
strain energy in a system free from geometric complications . The 
following set of nodal coordinates was selected: 

T 

{ui} = < Uj, Vj, Wj, Wi,g , Wj^g , Wi^sQ > i = 1,2, 3, 4 

(3.24) 

where s and 6 are the physical coordinates of any point on the element. 
Then using (3.1), (3.2), and (3.12) the transformation matrix [t] is 
constructed such that 

(ui) = [T] [aj . (3.25) 

The transformation matrix [t] is nonsingular and is easily 
inverted . 

Matrix [T] is given explicitly in Table 2 . 



G. STIFFNESS MATRIX 



The strain energy in the element is; 



U. = i 



[N] {c} r a ididri 

Substituting for stress [N] and strain [c] one gets: 
>T r^.l-,T^ rr r.,_,nT 



(3.26) 



Ug = i [ui)^ [E][P*]ro je a4Sn)[T'"] (ui) (3.27) 



where 

[P*] = [D] [P] 



(3.28) 
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The stiffness matrix is then: 

[K] = QJ [P*] [E] [P*]r o i 5^ St) [T"^ ] 



or 



[K] = [T-']'' [k] Cr"] 



(3.29) 



where [T]and [k] are both (24x24) matrices. 

Stiffness matrix [k] may contain some rigid body modes depending 
on problem type. However, by subjecting [K] to eigenvalue analysis, 
individual rigid body modes, if present, can be identified. 

When [k] does not contain all six components of rigid body motion, 
it must be modified accordingly. The method developed in Ref. 6 is 
applied as shown below: 



{U.} = [i ; R ] 



(3.30) 



where |u ^ j = < 6 , 6 , 6 , )S , )S , / Ll] is a (24x24) 



y z '^x y 

identity matrix, and [r] is the (24x6) matrix developed in Section D. 
Matrix [k] is then expanded to include all six rigid body modes; 

I . _ . . . . 

(3.31) 



[R]- 



[K] [: ; R] - 



where [R]"^ [K] [R] J is a (6x6) matrix. 

[k] : [k] [r] 
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ki’^ck] : cr:’'[k] [r] 
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(3.32) 



Solving for 
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(3.33) 



{u.^} = [[R]^[K][R] ] " [[k][r]] {uj 



However, if the stiffness matrix [k] already contains some of the rigid 

t T “1 

[R] [k] [R] will be singular. These 



modes have been previously identified and are merely suppressed in 
■|u^^ j- and in the corresponding column of [R] . Matrix [R]'^ [k] [R]J 



then becomes nonsingular. Solving (3.33) and substituting into the first 
(3.32) leads to the following expression: 

( 

{f } = (^ [k] - [[k][r]][Cr]^[k] [r]] ^[[R]^[k]]) { Uj } 

(3.34) 



The stiffness matrix modified to include rigid body modes is then: 



[k*] = [ K ] - [[K][R] ] [[R]"^ [k] [R]] " [[R]^[K]] 



Stiffness matrix [K*] is strain free for any rigid body motion. 
H. DEFINITION OF MATRICES 



In the following sections of this thesis, the stiffness matrix [K*] 
derived above is referred to as CSFE. The matrix [R]'^ [K][R] J 



developed in (3.31) is hereafter referred to as K22. For the general case 
where all six rigid body modes are absent, K22 is a (6x6) matrix. In the 
case of a cylindrical shell segment where four rigid body modes are 
absent, K22 is (4x4), and in the case of a conical segment where five 
rigid modes are absent, K22 is (5x5). 
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I. LOAD VECTORS 



Loads and moments acting on the element surface must be repre- 
sented by vectors which are consistent with displacements and rotations 
in the generalized coordinate system. 

A consistent load vector then, is defined as a vector made up of 
fictitious quantities which would give, after an inner product with the 
nodal point displacement, the same work as the real load. 

1 . Concentrated Load Vectors 



For concentrated loads of F 




of the element, the work done is 



W = F;.u + F V + F„w 

^ V C 



or more compactly 




Substituting for {uj from (3.12) and (3.24) gives 




The consistent concentrated load vector is then: 






LOAD 



(3.35) 



2 . Pressure Load Vectors 



The work done by a pressure load p is; 




A 



p w d A 



(3.36) 



then if 
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{q} = < 0 0 p > " 

direct substitution of known quantities into (3.36) gives: 

W=JJ{q3'^ {u} d^dT] = JJ{q}'^ pj rd4dT)^|ui| 

The consistent pressure load vector is then: 



{wl= t-’J'-dedt,. 



(3.37) 
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Table 1. Transformation Matrix [R] Between Nodal Coordinates 
and Polynomial Coefficients 
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where: 
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Table 2: Transformation Matrix LT] Between Nodal Coordinates 
and Polynomial Coefficients (Columns 1-12) 
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Table 2: (Continued) (Columns 13 - 24) 
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where : 
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IV. NUMERICAL CONSIDERATIONS AND PROGRAM INFORMATION 



A. GENERAL 

All numerical computations were accomplished on the IBM 360/67 
digital computer using FORTRAN 4 language. The entire program was 
coded in double precision arithmetic. An overlay scheme under control 
of the linkage editor was used to conserve core storage. 

B. QUADRATURE 

It was evident that the integrands of (3.27) and (3.37) could not be 
manipulated into closed form. Since the above expressions were in 
nondimensional form, lower and upper limits of integration zero and one 
respectively, it was advantageous to use Gaussian quadrature . An 
additional advantage of using Gaussian forms is that a polynomial of 
degree 2n-l can be exactly described by n points [9]. For polynomials 
of unknown degree, maximum accuracy is achieved by selecting a large 
number of points. Additional accuracy, however, is obtained at the 
expense of computation time. Investigation by the author showed that 
Gauss Legendre— 8 point quadrature provided accuracy within the range 
of double precision round off error for shells whose geometry could be 
described by analytic expressions. For shells of more complicated 
geometry, required accuracy is obtained when standard finite element 
mesh size reduction is accomplished. 
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C. LIMITING CASES 



In this thesis, cylindrical shells, truncated cones, and flat plates 
have been considered as special cases. When a shell is specified as 
one of the above forms, geometric approximation is precluded. 

Exact solutions for these shells provided a powerful tool for the 
analysis of the stiffness properties of the general case of the curvilinear 
shell finite element. 

D. CASE WHERE 0 IS EQUAL TO ZERO 

Singularities exist in the geometric approximation whenever the 
element is constructed in such a way that a point on its surface is 
locates at 0 equal to zero (see Figure 1). Discussion of the effects of 
this singularity is contained in Chapter 5, Section B, 2. 

E. INITIAL CONSTRUCTION OF THE STIFFNESS MATRIX 

The author originally constructed the stiffness matrix by constructing 
all matrices in (3.27) and then integrating the entire expression. 

Accurate results were achieved. However, the construction of the stiff- 
ness matrix for a single element required 500,000 bytes of computer stor- 
age and approximately 20 minutes of computation time. This was 
obviously unsatisfactory for other than pedagogical application. 

Fortunately, identical results were obtained by individual stiffness 
matrix element computation and by taking advantage of symmetry. This 
was done in the remarkably short time of 15 seconds using only 143,000 
bytes of computer storage. 
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CSFE was constructed and evaluated in a separate program prior to 



being incorporated in a direct stiffness program. 

F. DIRECTION STIFFNESS PROGRAM 

Numerical solution of engineering problems was accomplished 
through the use of a stress analysis program which was based on the 
direct stiffness method. The total stiffness matrix and the load vectors 
were assembled using programming techniques furnished by Professor G. 
Cantin. This program gave numerical values for displacements, moments, 
stresses, and strains at each nodal point of the mesh. A flow diagram of 
the computer program is presented in Figure 3. It can be seen that the 
program is functionally divided into two parts. In the CSFE portion, the 
stiffness matrix and its associated pressure load, strain-displacement, 
and stress-strain relationships are computed for individual elements. 
Results of these computations are stored on disk units. In the stress 
analysis portion, individual element information is retrieved from disk 
and assembled in a manner consistent with the direct stiffness method. 

The author's primary purpose for construction of the stress analysis 
portion of the program was to verify the stiffness properties of the 
curvilinear shell finite element. It was therefore, designed to handle 
problems of limited size. Specifically, the program is limited to a 
maximum of 36 elements, 49 nodal points, or 294 equations. Coding 
methods used, however, are general to allow for future program size 
expansion. The computer program is found in Appendix C. 
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V. DISCUSSION 



A. INTEGRITY OF THE CURVILINEAR SHELL FINITE ELEMENT 

Confidence in the integrity of CSFE was established by comparing" 

stiffness properties of limiting cases (flat plate and cylindrical shell) 
with properties found in the literature. Element by element comparison 
with values given by Bogner, et. al. [Ref. 10] revealed exact agreement 
for the case of the flate plate. In the case of the cylindrical shell 
element similar comparison was made with values given by Cantin Q^ef. 
l]; once again exact agreement existed. 

Confidence in the integrity of the stress analysis program was 
established by direct comparison of problem solutions using CSFE with 
solutions found in the literature. Results of comparison are discussed 
in the following sections. 

B. NUMERICAL SOLUTION OF CLASSICAL PROBLEMS 

The classical problems analyzed in this section were selected to 
demonstrate the integrity of the element, and to exhibit both the 
flexibility and limitations associated with CSFE. 

1 . Pinched Cylinder 

The pinched cylinder illustrated in Figure 4 has been analyzed 
by Bogner, et al [Ref. 11] using a (48x48) element and by Cantin 
O^ef. 6] using a (24x24) element. Table 3 is a comparison of results 
obtained using CSFE with those found in the literature. The obvious 
agreement of results confirms the integrity of CSFE. Moreover, it is 
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100 lbs 
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10. 3 5 i n. 
10500 ksi 
0.3125 



Fig. 4 PINCHED CYLINDER 
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evident that CSFE has the desirable characteristic of rapid convergence 



to the true solution (-0.1139) given by Cantin. 



No. of 
EQS. 


Bogner et 
(48x48) 


: al 


Cantin 

(28 X 28) 


CSFE 
(24 X 24) 




Mesh 


Displace . 

_Jin_.) _ 


Mesh 


Displace . 


Mesh 


Displace . 


54 






2x2 


-0.0931 


2x2 


-0.0932 


108 


2x2 


-0.0808 










150 






4x4 


-0.1126 


4x4 


- .1128 


180 


2x4 


-0.1098 










294 






6x6 


-0.1137 


6x6 


-0.1138 


486 






8x8 


-0.1139 






726 






10x10 


-0.1139 







Table 3. Displacement under load P for pinched cylinder 



2 . Spherical Cap 

The spherical cap illustrated in Figure 5 demonstrates the 
general case of CSFE. This problem which has an exact solution given 
by Timoshenko D^ef. 12] is useful for demonstrating the capabilities and 
limitations of the new element. Figure 6 contains a plot of meridional 
bending moments.^ 



^ The signs of Timoshenko's results have been changed to agree with 
the author's convention. 
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Fig. 5 SPHERICAL CAP 




Fig. 6 BENDING MOMENT M4- VERSUS-^ 
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It is worthy of note that the pole of the cap was located at a singular 
point (0 = 0). In order to fulfill connectivity requirements, a small rigid 
polar "cap" was assumed to exist. Boundary conditions were applied to 
allow this "cap" freedom of displacement in the z-direction at zero slope. 
Zero slope specification had the effect of introducing a false moment in 
elements adjacent to the polar cap. This false moment is visible in 
Figure 6. It is important to recognize, however, that the detrimental 
effect of the erroneous specification is localized in the near-pole 
elements . 

Results indicate that the bending predictions are quite accurate 
even though a coarse mesh was used. 

C. CONCLUSIONS 

A finite element for a shell of revolution which has arbitrarily curved 
meridians has been developed and tested. Element description was based 
on displacement functions which met all requirements to guarantee mono- 
tonic convergence to an exact solution by mesh size reduction. 

Results for classical problems obtained by using the new element 
indicated that good approximations for bending stresses could be obtained 
with a very coarse mesh. Mesh size reduction will be necessary to judge 
element performance in problems where membrane action is predominant. 

The geometrical deficiency where (0 = 0) is a limitation which must 
be recognized by users. Its effect is problem dependent. For the case 
of the spherical cap it was seen that the detrimental effect of artificially 
imposed boundary conditions at this singular point did not propagate 
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throughout the solution. On the other hand, the author was unable to 
successfully approximate stresses in a pressurized toroidal shell due to 
error propagated from similar artificial boundary conditions . Time 
precluded the determination of a method to efficiently handle this 
singularity. 

Not withstanding this geometric limitation, results indicate that the 
new element (1) is accurate, (2) requires low computational effort and 
(3) does not require an analytical description of shell geometry. Based 
on these facts it is concluded that the curvilinear finite element can be 
efficiently used in many engineering applications. 

D. RECOMMENDATIONS 
It is recommended that: 

> (1) the stress analysis program be expanded in order to handle 

realistic engineering problems; 

(2) a consistent gravity load vector be included in the stress 
analysis program; and 

(3) a study be conducted to determine a method of efficiently 
handling the singular point which exists where (0 = 0) in the new element. 
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APPENDIX A NUMERICAL EXAMPLE OF THE STIFFNESS MATRIX 



In this appendix, a typical example of the stiffness matrix CSFE 
is presented. Eigenvalues, before and after rigid body motion was 
included, and matrix K22 are also presented. Double precision figures 
were truncated to conserve space. 
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STl'^FNPS'; 'MATRIX FOR A TYDICAI ‘■.HFLL ELF^ENT 



YCUNG'S MOnULUS = 
PCI SSOK'S RATIC = 
RAOIU <^ 1 
RADIUS 2 

LENGTH = 

OHI 1 
PM I 2 

TMIC'-'N-SS 1 ' 
THIOXN-^SS 2 
A I, CM A - 



I .oooDon 

0.300000 
7.071100 
3 .G< 0300 
2.^13000 
0.7R5A00 
I .047200 
I .000000 
1.000000 
0.2M300 
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7.C173r-C2-3.77CCn- 


-01 5.44800-02 2.451 2P-0? 2 . 1 1 5 7 0-0 2-7 . n8 7 C-0 3 


1 .73^‘=^-C2-‘=.4^PC0- 


-C2-1 .01660-01-3.9 26 7 0-0 2-1 .0 7 69 0-01 3 .4^680-0 2 


1 .■^C*^5C-C2*-3 .451 20- 


-02-3.82679-02-1 .30450-02-3.70190-0? ^.81130-03 


2.6^6^C-C? Z.ll^^cp- 


-02 1.0'^990-ni 3.7on'^-02 8.19910-02-7.72570-02 


l.=f 1 ^r-C2-7.l9P70- 


-03-3.44680-02-6.811 3D-0 3-2. 72 5 70-0 2 2.53640-0^ 


1.99rfr>-ci 1.5^710- 


-01-3.73 770-0 2-1 .2 973 0-0 7-2.72 700-07 1 .5 3980-0 2 


l.-^p'^n-Cl-2.933^‘D- 


-01 5.14790-02 1.72290-02 7.30430-03-2.40^30-03 


?. '^9340-02-4 .4P12C- 


-02-1 .66 5 8 0-0 2-3.3 3340-0 2-4.4^^9 3 0-02 3.66670-0 2 


2.512*^0-02 1.8715D- 


-02 3.19620-02 4.12no-02 3 .89070-02-3 .03520-02 


3.^35^0-02 2.76?ad- 


-0? 3. 1 0660-02 3.23960-07 4 .70 p4p-o2- 3 . 1 57 20-02 


1.5*^*^7r-C2 2.5^^^D- 


-03 2.80870-02 2.40430-02 2.76510-02-1.68960-02 


(CCL. 7 TC 12 ) 
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- 1. 9 'JP60-C 1-1. 59710-01 -3. 73270-02-1 .2 8730-02 2.7 2700-0 2-1 .53980-0 2 



■ 1.9f» 8 70-01-2. 933AD-01 -5. 14900- 


-0 2-1. ■'2290-0? 7.304 30-03-2.494 30-0 3 


3.59’40-C2 4.48120-02-1.66530- 


-02-3.33340-02 4.44880-02-3.88670-02 


2. 51250-C2-1 .87150-02 3.19620- 


-02 4. 1215D-02-’ .89070-0? 3.03520-02 


3.83P60-C2 2.768AD-03-3.1036D- 


-02-3.23980-02 4 .79 240-02-3 . 1 57 20-0 2 


I. 65570-02 2.54940-03-2.80870- 


-02-2 .404 30-0 2 2.76 510-02-1 .68880-02 


4.73250-02 7.C1730-02 1.73650- 


-07 1.70550-02-2.68660-02 1.56130-02 


7.01730-02-3.77000-01-5.44800- 


-02-2 .451 20-0 2 2.11 59 0-0.7-7.19870-07 


1.73650-02 5.44600-02-1.01660- 


-01-3.82670-02 1.07890-01-3.44680-02 


1.70550-02 2 .45120-02-3. 326-'0- 


-02-1 .30450-02 3.70190-02-6.81130-03 


2. 68660-02 2.11590-02-1 .07fl«0- 


-01-3. 701 <5 0-07 8.19 810-02-2.77 570-0 2 


1. 56130-02-7. lce7D-03 3.-^4680- 


-02 6. 8113D-03-2. 72570-02 2.53640-03 


3. 45900-01 1.94390-01 5.43790- 


-02 2.60770-07-3.49390-02 1.86310-02 


1.94390-01 4.77170-01 7.45310- 


-02 2.38510-02-2.29220-02 7.56010-03 


'^.4^790-02 7.48310-02 1.86280- 


•01 1.33400-01-1.54730-01 7.42140-0? 


2.60770-02 2.38510-02 1.33400- 


•01 1.7145D-01-1 .03860-01 7.57310-02 


3. 4 9 390- 0 2- 2. 2 <=2 20-02- 1.54730- 


•01-1 .03860-01 2.01990-01-8.84570-0? 


1.86310-02 7.56010-03 7.42140- 


•02 7. 57310-02-8. 84570-02 a. 67230-02 


1. 53090-01- 5. 422 10-03-1. 85 260- 


•02-2.36640-02 3.04170-02-1 .73980-0? 


4.67<0D-02 1.70660-01 2.69420- 


•02 1.80400-02-5.14550-03 2. 05150-03 


4.16120-02 2.00610-02-2.98750- 


•02-4.40050-02-1.79460-02 8.92330-03 


2.39650-02-1 .27510-02 3.75500- 


•02 2.68950-02 7.06530-03-5.55010-03 


3.811 2D-02-6.56C5D-04-1 .01370- 


•02-1.06240-03 7.0877D-02-1 .88300-0? 


1.73170-02-1.97280-03-5.40180- 


■03-3.31520-03 2.10610-02-1.30960-02 


(COL. 13 7C 18 ) 
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-2.f-C°CC-C? 2.9f^f'90-02 5.96220-03 1 .^0020-02 3.43610-02 1.54690-02 



2.<^66<;r-c2-‘^.03oin- 


-01-A .24970-02 1.S9000-02 3.74840-03 1. 29580-01 


5. =^2;n-C3 M*24<^cn. 


-02-1.632^0-01 5.01910-02 1.31320-01 3.96740-02 


1 .'IC02C-C2-1 .5Q0CP- 


-02 01 Q7q-o 2-2 .508 30-0 2-^ .35 740-0 2-1 .26150-0? 


3,^0MP-C2 3.^^04O* 


-C3-1 .31320-01 4.357^0-02 7.09810-02 2.959C0-02 


l.^^ecc-C2 1.2P‘=^8D- 


-03-3. °6 *^40-02 1.26150-02 2 . ’^590^-02 5.84220-07 


1 . coo^p-ci 1 .^p07p. 


-01 3.59340-02-2.51250-02-3. 80-02-1 .65570-0? 


1 . ^ '7 1P-C1-2.<5?34C- 


.Ol-^*.i^0l2n-O2 1.671^0-02 7.7^pAr-o3 2.54q4D-07 


?.7^?in-C2 


-02-1 .66580-02 3.1^620-02 1.10860-02 2 .8087P-02 


1.2^73r-C2 1 .72?^P- 


-02-3. 33340-02 4.12150-02 3.23^30-02 2.40430-0? 


2.72 7CO-C2 "^.3CA3P- 


-03-4.448R0-02 3.39070-02 4.7924C-02 2.7651C-0? 


l.f 3^PC-C2-2.47-^ 3P- 


-03 7 .P867Q-.02-3 .015 20-0 2-1 .1572^-02-1 .68 8PP-02 


i.*=3c<;n-ci ^.67*vcp- 


-02 4.161 2D-0 2-2 .396 5 0-0 2-7 .a 1 120-02-1 .73170-02 


5.^-221P-C3 1 .7C640- 


-O’ ?. 00 610-0 2-1 .2751 0-02-6.56050-04-1 .97280-07 

r 


l.P526n-C2 2.6^420- 


-C2 -2. 98 75 0-0 2 3 . 75 5 0 0-0 2-1 .0 1 370-02-5 . 40 1 03 


2.366 ^C-C2 1 .BO^CD- 


-C2-4.4 0C5D-0 2 2 .6 89*“ D-0 i-l .06240-0 3-3 . 3 15 20-Q7 


3.C^l70-C2-5.l^550- 


-01-1.79460-02 7.06530-01 2.087?c-02 2. 10610-02 


K”^3^?C-C2 2.C51PD- 


-03 8.92 33 0-03-5.5 501 0-0 3-1 . 8« 30C-0 2-1 .80960-0 2 


^.B6l lC-Cl-1 .5-29P- 


-01-6. 86000-07 2.^4990-02 4. 04830-02 1.7844P-0? 


1.^92^C-Cl 5.34C7D- 


-Cl 7.28130-02-2 .72040-02-^.53560-03-1.99200-03 


6,8^CCC-C2 7.2^130- 


-02 2.46880-01-1 . 38200-01-1 .7 ?o i q-o 1- 7 . 1 3 1 70-0 ? 
1 


2.9^cc;[:-C2-2.3204n- 


-02-1.38200-01 1.626P0-O1 8.09860-02 6.10670-02 


4.0'^fl;af>-C2-6.B356n- 


-C3-1 .72910-01 8.89360-02 1.7341C-01 6.46350-02 


U7.qA^C-C2-1.992CP- 


-03-7.13170-02 6.1067D-02 6.44350-0? 4.6593C-0? 


(CPL. 19 TC 24 ) 
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£IGF-NVALl = S flPFTRF PIGIC RCDY ''CTICN INCLUOEC 

1.6C7'?r CC (?. 735^0-01 8.03»'^0-01 7.50900-01 5.92120-01 '.9129C-01 
5.57^Cn-Cl 5.12150-01 1.69550-01 1.01500-01 9.7fl37C-0'> 8.5<?P50-02 
^.I0ccn_c2 5.93«‘'<O-02 1.972'»0-02 1.550^<C-0? 1.7?n00-02 6.22710-03 
5. '^1600-03 3.650^0-0'^ ?. 62050-05 3.50450-05 1.00450-06 1.55470-07 



VATRIX K22 

2.15] 20- C 6-5.54 ccr -O':-! .54 4B0-04 5.5 53 5 0-0*: 1 .50 77 0-06-2 .0'5 30C-0 5 
6.5 45 50-C5 5.621 5C-05-2.172in-0 5-6.5 7&O0-0 5-5.56 35C-05 1 .*:oo5C-0 3 
1 .66 5 CO-C6-3.1721D-C5 2.1595 0-05-1 .0^6 2 0-0 ‘t 9 . 2 5Q9 r-0--6 . 505 2C- 1 7 
3.56 3 6D-C5-6.5 06SI?-06-! . 085? 0-05 5 . 3 25 3 0-0 3-3 . 5*’ 9 3C-05-5 . 1 5 250-0 3 
1.657 70-06-6.55360-05 9. 25 090-05-9 . 6 03 3 0-0 v 1 .23170-02-9.11300-06 
2 .0-:: 3Cf'-C5 1 . 55 C5r-0‘--2 . 63 5 90-1 6-6 . 1 526 0-0 3-B . 1 1 30f'-07 P . =0 7 90-0':' 



=lGcNV3lLF£ ftPTFR RIGID POOY MOTION WAS INGLUCEO 

1.6C3Sn CC 8.7532C-01 7.03650-01 5.75820-01 5. 70550-01 5.6P020-01 
3.765 S0-C1 7.B611D-01 1 .03530-01 l.OOS^O-Ol 9.77770-02 7.2P'»10-07 
5.22770-C2 7.06990-02 1.75250-02 1.50520-02 1.17500-02 32370-03 
F.2':'71C-n 6.11250-17 3.31350-1 '’- 7 . 7 92 3 0-3 a - 4 .55 1 sc -1 7 - 7 . 7 30 10- 1 7 
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APPENDIX B COMPUTER OUTPUT 



In this appendix, a complete set of results for the spherical cap 
analyzed in this study are presented. Results for the spherical cap were 
obtained using a 4x4 mesh. The mesh layout accompanies the problem. 
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1 

2 

2 

42 

0 




SPHFRICAL SH6LL WITH EXTERNAL PRESSURE (1 PSI) 
TOTAL NUMBER ELEMENTS 

NUMBER OF OIFFEt’PNT TYPE ELEMENTS = 

TOTAL NUMBER CF NODES 

NUMBER nc NOnSS WITH BGUNOAPY CONDITIONS = 
HALFBAN’D WIDTH = 

NUMBER CF POINTS WITH CONC ENTP 4T EO LOADS = 
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FL 


SPEC 


F /nnis 


PAOl/PAP? 


PHI 1/PHI2 


HI /H2 






L /ALPHA 


1 


GEN 


4.00000 OfS 


7.14^10 


00 


3,'4^100-0? 


■’•OOOCD 


00 


1 


.20470 01 






1 .(S6700-01 


1.562^0 


01 


1.74530-01 


3,00000 


00 


3 


,^2^cn-ci 


2 


G«^N 


4.00000 04 


T.562S0 


0! 


1 .7 70-0' 


3,000Cr> 


CO 


1 


.33410 n 






I •66700-01 


?,3*"57D 


01 


3,22090-01 


3,00000 


00 


3 


.927CO-C1 


"3 


GEN 


4.COCOO 06 




01 


3,22«9n-01 


^•OOOCD 


oc 


1 


, 335 1^^ Cl 






1 •fb-’OD-Ol 


A .OPSoo 


0^ 


A .-^1240-01 


”^,00000 


00 


7 


,G27C0-C1 


1 


GPN 


' ,00000 06 


A .08500 


01 


' •“'12^0-0! 


^ .00000 


00 


1 


.2467^ 01 






! ,6^ TOO-01 


,16 2 70 


01 


6 • 10870-01 


^ .00000 


00 


■3 


,927Cr-C^ 



CC^^ECTI V! TV 
T 


J 


V 


1 


Type 


1 


7 


7 


4 


1 


/ 


7 


>) 




9 


7 


4 


1 


0 


7 


L 


S 


' 0 


Q 


/ 


6 




12 


1 ! 


1 


7 


0 


13 


1 7. 


9 


P 


Q 


14 


13 


Q 


c 


10 


IS 


1 A 


f 


11 


12 


1 7 


1 4 


1 


12 


1 3 


1 8 


\ ^ 




1 ^ 


14 


1 o 


1 r 


7 


1 A 


15 


7.0 


1 c 




1 P 


1 


22 


21 


1 


17 


14 


23 


22 


9 


1 8 




24 


7 3 


7 


1C 


20 


2 4 


74 
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BOLNDAPY CONDITION'; 



JOINT 


X TRANS 


Y TRANS 


Z TRANS 


W, X 


1 


0 


0 


1 


0 


A 


0 


0 


1 


0 


11 


0 


0 


1 


0 


16 


0 


0 


1 


0 


21 


0 


0 


1 


0 


5 


0 


0 


0 


0 


10 


0 


0 


0 


0 


15 


0 


0 


0 


0 


20 


0 


0 


0 


0 


2 = 


0 


n 


0 


0 


2 


1 


0 


1 


\ 


•3 


1 


0 


1 


1 


4 


1 


0 


1 


1 


'>? 


1 


0 


1 


1 


2? 


1 


0 


1 


1 


24 


1 


0 


1 


1 


7 


1 


1 


1 


1 


P 


1 


1 


1 


1 


C 


1 


1 


1 


1 


12 


1 


1 


1 


1 


18 


1 


1 


1 


1 


14 


1 


1 


1 


1 


17 


1 


1 


1 


1 


1? 


1 


1 


\ 


1 


1^ 

0 1 $ P L 6 r 


1 


1 


1 


1 



W,Y 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

0 

0 

0 

0 

0 

0 

c 

0 

0 

c 

c 

0 

0 

n 



0 



)<Y 

0 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

c 

c 

c 

c 

c 



jriNT L V 

1 -c.c -o.c 

Z ’.C^UO-C5 -o.c 

^ (".4^3CO-0^ -0.0 

^ f . ic'H^n-c*' -0.0 
-c.c -0.0 

^ -C.O -0.0 

7 3.C-120-05 -3.?0‘5in-17 

« f.^33CP-CS 2 
c; (^. l«i'>4n-05 l.TSqcrj-iA 

10 -c.c -0.0 

11 -c.c -0.0 

12 3.C^120-C= -A.^^ftOD-l-' 

13 6.40300-05 ^.30530-1*^ 

1^ 6.15240-0= 2.2'3^jPo-is 

1= -c.c -0.0 

I-'- -c.c -0.0 

17 3.Cn20-05 -3.20740-17 

IB 6.4?3C0-0= 1.99380-14 

1C 6. 15240-05 1.75930-14 

2C -C.C -O.C 

21 -C.C -0.0 

22 ?. 09120-0= -0.0 

23 6.^83C0-0= -0.0 

24 4 . 1 ^ 240-05 - 0.0 

2= -C.C -0.0 



W W,X W,Y 

-4.00210-04 -0.0 -O.C 

-5.6^1 10-04 -5.1250f'-06 -C.C 
-4.3'>31O-04 -1.50350-05 -0,0 
-1.82260-04 -?.15®3'"-0= -0.0 
- 0.0 - 0.0 - 0.0 
-5. 00210-04 -0.0 -O.C 

-*^.64110-04 -5.12 = 00-04 -0.0 
-4.323n-04 -1.50350-05 -0.0 
-1.32240-0^ -2.15^30-05 -0.0 

- 0.0 - 0.0 - 0.0 
-4.00210-04 -0.0 -0.0 

-5.6^1in-04 -5.12500-06 -0.0 
-4.32310-04 -1.50350-0= -C.O 
-1.8224D-04 -2.1=8^0-05 -0.0 
- 0.0 - 0.0 - 0.0 
-4. 00210-04 -0.0 -0.0 

-5.44110-04 -5.12500-04 -0.0 
-4.32310-04 -1. 50350-05 -0.0 
-1.82260-04 -2.19830-0= -0.0 
- 0.0 - 0.0 - 0.0 
-4.00210-04 -0.0 -0.0 

-5.44110-04 -5.12500-06 -0.0 
-4.32310-04 -1.503=0-05 -0.0 
-1.82260-04 -2.15830-05 -0.0 
- 0.0 - 0.0 - 0.0 



vj , XY 

- 0.0 

- 0.0 

- 0.0 

-o.c 

- 0.0 

- 0.0 

-C.O 

- 0.0 

- 0.0 

-o.c 

- 0.0 

-o.c 

-C.C 

-o.c 

-C.C 

-o.c 

-o.c 

- 0.0 

- 0.0 

- 0.0 

-o.c 

- 0.0 

-C.O 

-c.c 

0.0 



55 



STRFS5F? 



JCINT 
1 -1 
2 -A 

^ -1 
7 

P -c 
C - c 
10 -6 
11 -1 

12 -A 
1? 

14 -c 
1*5 -6 
16 -1 
17 -A 

JQ -C 
1C -5 
20 -6 
21 -1 

22 -A 

23 

->L 

-6 

T^ < 

JCINT 

1 

2 -3 

•> _ ■a 

/ -4 , 

c .4 
4 - C ^ 

7 

Q - 3 , 
G -4 , 
1C -4, 
1 ! > 
12 

13 

-A, 
IF -A, 
-G, 
17 

IP -3, 
IQ -A, 

20 -4. 

21 -c. 

22 -3, 
2^ -3, 

'>L -4 , 

- 4 , 



KX KV NXV »^X wy mxy 



,32'^‘;D 02 -1.0200D 02 
,';A17D 01 -6.00750 01 
.0A3AD 01 -A.0216D 01 
,P23CD 01 -1.790BD 01 
.CA29D 01 -1.0073D 01 
,32^50 02 -I .02000 02 
.GA17D 01 -b. 00750 01 
.CA3AD 01 -A. 02160 01 
.P2300 01 -1.7Q0P0 01 
.04290 01 -1.00730 01 
.32990 02 -1.02000 02 
.9A1 70 01 -6.00''''0 01 
.CA240 01 -4.02160 01 
.°23CO 01 -1.79030 01 
.06290 01 -1.00730 01 
. '^’GGO 02 -1 .02000 0? 
.CA170 01 -6.00750 01 
.0^*340 Cl -4.0216D 01 
. P2TCD 01 -1 .790P0 01 
.04290 01 -1.00730 01 
.37990 02 -1 .02000 0^ 
.CA1*»D 01 -6.00750 O! 
.04340 OT -4,02160 01 
.«23CO 01 -1.7Q0B0 01 
, CA2QO 01 -1 ,00730 01 



SNX SKY 

66 = PD-C6 -6, 6*^250-06 
2S3fD-06 -A.31QP0-06 
>A42P-05 -2.^ = 070-06 
6C370-06 -6.53440-07 
P95P0-C6 2.99520-22 

66560-0*- -6,6 = 250-06 
25360-0= -4.31950-06 
6AA20-06 -2.65070-06 
60770-06 -6.53440-07 
pc«;Qn-06 -3.2250O-2‘> 
=6560-06 -6,65250-06 
25360-06 -4,31950-06 
=4420-04 -2.65070-06 
60370-06 -6. 93440-07 
P95P0-06 -3.32500-22 
66550-06 -6.65250-06 
25360-06 -4.31950-06 
6442D-06 -2.65070-06 
60370-06 -6. 93440-07 
P95PO-06 -3.32500-32 
66560-06 -6,6=250-06 
29360-0= -4.31950-06 
64420-06 -2.65070-06 
6C37D-06 -6,83440-07 
59=60-06 -9.6=130-27 



0.0 -4.65920 00 

-3,50500-12 -4.44650 00 
-3.56360-11 -7,62730 00 
-5,10700-11 -2.6T6S0 00 
0.0 3.17990 01 

-1.311 60-11 -4.65920 on 
-1.37460-1 1 -4. '+4650 00 
7.35530-12 -7.62730 OO 
-4.47250-11 -2.61650 00 
-7,20040-1 1 3 , 1 7Q90 01 

-1.97090-1' -4,65370 00 
-2.11660-11 -4.^4650 00 
7.00610-11 -7.627?0 00 
-3. 07170-1 1 -2.61650 00 
-9,39660-1 1 ^1 

-1.31260-11 -4,65920 CO 

-1 .37630-1 1 -4.44*<^o 00 
7,32590-12 -7.627^0 00 
-4,474^0-11 -2.61650 00 
-7.19'570-n 3.17990 01 

0,0 -4,65920 00 

-3,91210-12 -4.44680 00 
-3.87030-11 -7.62730 00 
-5.11130-11 -2.616=0 00 
0.0 7.17990 01 



SKXY S^x 

0,0 -=.16=50-07 

-7.39670- 1<3 -4.29550-07 
-”.51320-15 -7. '+4^70-07 
-9.92040-18 -2.06650-07 
0.0 3 . -*35 10-0 = 

-2.5=0=0-15 -=.06=50-0^ 
-2.67300-U, -4,^9050-07 
1 ,43010-1 5 -7.^*4670-07 
-5.69740-lP -2.16680-07 
-1.40010-17 3.47510-0= 

-3.53240-18 -5.1655P-07 
-4.1 1 590-1 '' -4,230=0-07 
3.90040-15 -7.4-*67D-07 
-5,07310-15 -2.066P0-0” 
-1 .52720-17 7.-^3510-06 

-2.5=230-15 -=,16=50-07 

-2 . 67670-15 -4.295 = 0-07 
1.42390-18 -7.44670-07 
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APPENDIX C COMPUTER PROGRAM 



In this appendix the stress-strain analysis program is presented. 
User and functional information are presented through the use of comment 
cards located throughout the program . 

The program input data deck is illustrated below. Data values 
presented were those used for the spherical cap (2x2 mesh) . 
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C STRESS ANALYSIS PROGRAM FOR CONSTRUCTING ANP TESTING CSFE 
C 

C CONTROL SUBPROGRAM 

IMPLICIT REAL*8 (A-H,0-Z) 

,RFAL*8 K1,K2,KT,L 

COMMON/Bl/‘=,Hl,H2, POI S,RAD1 ,RAD2»Lt PHI 1 » PHI? , ALPHA, 
«THETA,KASE 

C0MP0N/MBN/8L0AD(29A) ,NRC(50,7J ,NCON(36,5) 
COMMON/PARM/NEL ,NELT,NDPT, NPBC, NBAND.MFO 
50 DO 100 1=1,294 
100 3L0A0( I )=C.ODO 
CALL INPUT 
IFINELT.EC.O ) STOP 
CALL OVl 
CALL 0V2 
GO TO 50 
ENP 



SUBROUTINE OVl 

IMPLICIT REAL*a (A-H,0-Z) 

common /MRN/BLOAO( 294 ), NBC ( , 7 I , NC ON ( 3 6 , 5 ) 

COMMON/PARM/NEL ,NELT,NDPT, NPBC, NBANO,NFQ 
DIMENSION T0TK(294,54) 

DO 100 1=1,294 
DO 100 J=l,54 
100 TOTK( I,J)=C.CDO 
CALL FORML 
CALI MERGER (TOTK) 

CALL BANSOLI TOTK, B LOAD, NEO,NB AND) 

RETURN 

END 



SUBROUTINE 0V2 
implicit REAL*8 (A-H,n>Z) 
CALL OUTPUT 
CALL PSSN 

return 

END 



SUBROUTINE INPUT 
IMPLICIT REALMS (A-H,0-Z» 

REAL*8 K1,K2,K?,L 
INTEGFP*4 SPTY(4),SPEC 

DIMENSION STFM0D(24, 24) ,T INV( 24, 24) , TINVT( 24, 24) 
DIMENSION PVECI24) ,CF0RCE(6) 

CQMMON/Pl/E, HI, H2, PO I S , RAD 1 , R AD 2 , L , PHI 1 , PHI2 , ALPHA, 
«THETA,KASE 

C0MM0N/MBN/BL0AD(294) ,NBC(50,7) ,NC0N(36,5) 
COMMON/PARM/NEL ,NELT,NDPT,NPBC, NBAND,NEO 
COMMON/SS SN/SS (24,24) ,SN(24, 24) 

DATA SPTY/'GEN CYL PLT CONE’/ 

C BLGAD IS TOTAL LOAD VECTOR 

C NBC IDENTIFIES THE NODE POINT WITH ITS BOUNPAPY CONDITION 
C NCON IDENTIFIES ELEMENT NODAL POINT AND ELEMENT TYPE 
C NRAND IS FCUAL TO HALFBAND WIDTH 

C NPCL IS THE NUMBER OF POINTS WITH CONCENTRATED LOADS 
NTK = 0 

C NTK IS CYLINDER TRACK COUNTER 
WRITE(6,1100) 

1100 FORMAT! IHl,//) 

READ (5,1000) 

WRITE (6,1000) 

1000 FORMAT(ftOH 

H ) 

READ (5,2000) NELT , NEL , NDP T , NP BC , NB AND , NPC L 
2000 FORMAT(6I10) 

IF(NELT.EC.O) RETURN 
NEQ=6*NDPT 

C NELT IS TOTAL NUMBER OF ELEMENTS 
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C NEL IS NUf^PER OF DIFFERENT TYPE ELEMENTS 
C NDPT IS TOTAL NUMRFR OF NODES 

C NPRC IS number of NODE POINTS WITH ROUNDAPY CONDITIONS 
WRITF (6,2500) NELT,NEL,NDPT,NPRC,NBANO,NPCL 
2500 FORMAT! /, 3BX, • TOTAL NUMBER OF E LEMF NTS • , 1 R X , • , 1 1 0 , / , 
«38X, 'NUMBER CF DIFFERPMT TYPE E L EME NTS * , PX , • = • , 1 1 0 , / , 
#3RX, 'TOTAL NUMBER OF NODES ' , 2 1 X , ' = ' , 1 1 0 , / , 

«3BX, 'NUMBER OF NODES WITH BOUNDARY CONDITIONS =',110, 
«/,38X, 'HAL PR AND WinTH',28X,' = ', I10,/,38X, 

W'NUMREP CF POINTS WITH CONCENTRATED LOADS =',I10) 
WRITF (6,3000) 

ROOO FORMAT (IHl ,///, 32X, ' EL SPEC E /POIS RADI/', 

#'RAD2 PHI)/PHI2 HI/H2 L /ALPHA') 

DO 100 1=1, NPL 

READ(8,R5nO) RADI, RA 02, L,PHI1,PHI2, ALPHA 
35CO FORMAT! 7F1C. 0) 

READ (5,AOOO) F , PO I S , H I , H2 , P , K A SC 
4000 F0RMAT(5F10.0, I 10) 

P IS THE PRPSSUPE LOAD 
KASE SPECIFICATIONS FOLLOW 

KASF=0 CURVILINEAR SHELL 
KASF=1 CYLINDRICAL SHELL 
KASE=2 FLAT PLATF 
KASE=3 TRUNCATED CONE 

IF A FLAT PLATE IS SPECIFIED WIDTH IS WRITTEN FOR ALPHA 
I F( KASE. EQ.2 ) AL PH A= 1 . 0 DO/ THET A 
IF(KASF.FO.O) SPEC=SPTY(1) 

IF( KASE.pq. 1 ) SPEC=SPTY(2) 

IF(KASr.FQ.2 ) SPPC=SPTY(3) 

IF! KASE. FQ. A) SPEC=SPTY(4) 

WRITE (6,4100) I,SPEC, E,RAD1,PHI1,H1 ,L ,PPI S,RAD2,PHI2, 
«H2, ALPHA 

4100 FORMAT !/ / ,32X, I 2,3X,A4,2X,1P5D12.4//43X,1P5D12.4) 

CALL ELSTF(STFMOD) 

NTK=NTK+1 

CALL WRDI SK(NTK,STFMnO,576) 

CALL TFCRM(T INV, TINVT ) 

CALL STRFSS(tinV,SS, SN) 

NTK=NTK+1 

CALL WRriSK!NTK,SN,576 ) 

NTK=NTK+1 

CALL WRnSK(NTK,SS,576) 

NTK=NTK+1 

CALL PLCADIP ,PVFC) 

CALL WRriSK(NTK,PVEC,24) 

100 CONTINUE 

WRITE (6,4400) 

4400 FORMAT (1 HI, /// ,32X, 'CONNECTIVITY') 

WRITF (6,4S00) 

4 50C FORMAT! //, 33X ,' I ',15X,'J',15X,'K',18X,'L',13X,'TYPE') 
DO 200 I=1,NELT 

READ (5,8000) ( NCON ( I , J ) , J= 1 , 5 ) 

8000 F0RMAT(5I10) 

WRITF(6,5500) ( NCON ( I , J ) , J = 1 , 5 ) 

5500 FORMAT! /,32X,I2,14X,I2,14X,I2,14X,I7,14X,I2) 

200 CONTINUE 

WRITF (6,6000) 

60C0 FORMAT! IHl ,/// ,?7X ,' ROUNDARY CONDITIONS') 

WRITE! 6, 6100 ) 

61C0 FORMAT!// ,27X, ' JOINT ', 3X, 'X TRANS', 5X,'Y TRANS', 5X, 

#'Z TRANS' ,8X,'W,X' ,0X,'W,Y',9X,'W,XY') 

C READ BOUNDARY CONDITIONS 
DO 300 I=1,NPBC 

READ (5,65CO) I NBC ( I , J ) , J= 1 , 7 ) 

6500 F0RMAT(7I10) 

WRITF!6,7000) ( NBC ( I , J ) , J= 1 , 7 ) 

7000 FORMAT!/, 2RX, 12, 8 X , I 1 , 1 1 X , I 1 , 1 1 X , I 1 , 1 1 X , I 1 , 1 1 X , 1 1 , 
#11X, II ) 

300 CONTINUE 

IF(NPCL.EQ.O) RETURN 
WRITE (6,T050) 

7050 FORMAT! IHl ,/// ,27X , 'CONCENTRATED LOADS') 
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WRITE (6,7100) 

7100 FORMAT (// ,27X, • JOINT *3X, 'FORCE X',5X,'FnRCF Y',5X, 
# 'FORCE 7',RX,'W,X',RX,'W,Y',9X,'W,XY') 

C READ CONCENTRATED LOADS 
00 500 I=1,NPCL 

READ (5,72C0) JNR , (CFORCF ( K ) ,K = 1 ,6 ) 

7?00 FORMAT (I10,6F10,0) 

WRITF (6,7300) JNR , ( CFORCE ( M ) , M=1 , 6 ) 

7300 FORMAT) /, 19X,I 12, 1P6D12.4) 

IB=6*JNR-6 
on AOO J=l,6 
IRL=IB+J 

400 8L0AD( IPL)=CFORCE( J) 

500 CONTINUE 
RETURN 
END 



FUNCTICN COSPHI (X) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*« K1,K2,K3,L 

C0MM0N/R1/E,H1,H2, POT S,RAD1 ,RAO?,L ,PHI 1 , PH I 2 , ALPHA, 
#THETA,KASE 

C0MM0N/B2/CPHI 1 , CPHI 2 , CO , C 1 , C2 
COSPHI =C0+C1*X+C2*X**2 

IF(KASE.GE.l ) C0SPHI = 0.000 
IFIKASE.FQ. 3) COSPHI =OCOS ( PHI 1 ) 

IF(DABS(COSPHI ).GT.1.0DO) GO TO 1 
RETURN 
1 CONTINUE 

CALL FLAG(2,ei00,C200) 

100 CONTINUE 
200 CONTINUE 

C VALUE ASSIGNED IS INTERIM VALUE ONLY 
C0SPHI=1 .000 
KASE=9 
RETURN 
END 



FUNCTION SINPHI ( X) 

IMPLICIT PEAL*8 (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

S INPHI = nSORT( 1. 0D0-( COSPHI ( X) )'<‘*2) 

RETURN 

END 



FUNCTION R(X) 

IMPLICIT REAL*8 (A-H,0*7.) 

RFAL*8 K1,K2,K3,L 

COMMON/ Bl/E, HI, H2, PO I S , R AD 1 , R AD 2 ,L , PHI 1 , PH I 2 , AL PHA , 
«THETA,KASE 

C0MM0N/B2/CPHI 1, CPHI 2, CO, Cl, C2 
IF(KASE.E0.2 ) GO TO 1 

R =RAD1 + L*(C0*X+(C1*X’«'*2 )/2 . 000 + ( C 2 *X* *3 ) / 3 ♦ ODO ) 

IF(KASE,EQ. 1 ) R=RAD1 
RETURN 
1 CONTINUE 

C VALUE ASSIGNED IS INTERIM VALUE ONLY 
R=1 .ODO 
RETURN 
END 



FUNCTION C(X) 

implicit REAL*8 (A-H,0-7) 

REAL*8 K1,K2,K3,L 

CCMM0N/B1/E,H1 ,H2, P0IS,RA01, RAD2,L, PHI1,PHI2,ALPHA, 
#THETA ,KASF 
Q=1.0D0/R (X ) 

IF( KASE.EQ.2 ) Q=0.0D0 
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RETURN 

END 



FUNCTITN OC(X) 

implicit RE/^L^e (A-H,n-Z) 

REAL’*'R K1 ,K2,K3,L 

COMMON/Pl / '=,H1 ♦H2,P0IS,RAD1 ,RA07,L, PHI I , PH I ? , AL PHA, 
«THFTA,KASF 
QO = C( X > /At PHA 
IF(KASF.FQ,2) OQ=THETA 
RETURN 
FND 



FUNCTION Rl( X) 

implicit REAL*8 (A-H,n-Z) 

REAL*R K1,K2,K3,L 

COMMON/Rl/E,Hl,H2t POIS,RADl,RAn2,L ,PHI 1 , PH I 2 , AL PH A , 
#THFTA,KASF 

COMMCN/B2/CPH1 1 ,CPHI2 ,COtCl ,C2 
IFIKASF.GE.l ) GO TO 1 
R1=-L*S INPHI (X ) / (C1+2.0D0*X«C2 I 
C THE FOLLOWING CARD PROTECTS AGAINST OVERFLOW IF X IS 
C EXACTLY AT A POINT OF INFLECTION 
IF(R1 ,GT. 1 .0C50 ) GO TO 1 
RETURN 
1 CONTINUE 

C VALUE ASSIGNED IS INTERIM VALUE ONLY 
Ri=i.ono 
RETURN 
END 



FUNCTION Q1(X> 

IMPLICIT RFAL*3 (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

COMMON/Rl/P,Hl,H2,POI S,RAD1 ,RAD2,L,PHI 1 ,PHI2 , ALPHA, 
<#THETA,KASE 
Q1=1.0DC/R1 ( X> 

IF(KASF.GE.l) Ql=O.ODO 

RETURN 

END 



FUNCTION ORIX) 

IMPLICIT REAL+B (A-H,0-Z) 

REAL*8 Kl,K2,K3,L 

COMMON/P1/e,h 1,H2, POIS,RADl,RAD2,L, PHI 1 ,phT2,ALPHA, 
#THETA , KASE 
DR = L*COSPHI I X) 

RETURN 

END 



FUNCTION DRKX) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

C0MM0N/B1/E,H1,H2, POISfRAOl ,RAD2,L , PHI 1 , PHI? , ALPHA, 
#THETA, KASE 

COMMON/B2/CPHI 1 , CPH I 2 , C 0 ,C 1 , C2 
IFIKASE.GE.I ) GO TO 1 

nRl = DR(X) / SINPHI (X)+2.0D0»C2*(-R1(X ) ) / I Cl + 2. OQO+X -^C? ) 
RETURN 
1 CONTINUE 
ORl=O.CDC 
RF^URN 
END 



61 



SUBROUTINE ELSTFI STFMOD ) 

C ELSTF IS THE CONTROL PROGRAM FOR CONSTRUCTING THE 
C FLEMENT STIFFNESS MATRIX. 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

DIMENSION TINVTI 24,24) , T I NV { 24 , 24) , XG { 8 ) , W T W ( 8 , 8 ) , 
#STIFF1 (24,24),STF( 24,24),STK(24, 24) ,WT( 8), 
#STIFFF(24,24) , STF MOD ( 24 , 24 ) 

COMMON/P 1/E, HI, H2, PO I S, RADI ,RAD2 ,L ,PHI 1 , PHI2 , ALPHA, 
#THETA,KASE 

COMMCN/B2/CPHI 1 , CPH 1 2 , CO , Cl , C2 
CPHI1=DC0S(PHI1 ) 

CPHI?=DCCS(PHI2 ) 

CO=CPHI 1 

Cl = 2.0nC*( 3.CD0‘>=(RAD2-RA01 )/L - 2. OOO^C PH 1 1 -C PH I 2 ) 

C2 = 3.0D0+(-2.0D0*( PAD2-R AD 1 ) /L + CPH 1 1 +C PH 1 2 ) 

IF( KASE.EQ.O) GO TO 20 

ci=o.ono 

C2=o.ono 

20 CONTINUE 

DO 50 1=1,24 
no 50 J=l,24 
50 STK ( I , J )=0.0D0 

C ZEROS OF GAUSSIAN 8 POINT FORMULAS ARE 
XX 1 = . 960 288 8 564975362300 
XX2=. 79666647741 36 2674D0 
XX3=. 5255324099163289900 
XX4=, 1834346424956498000 
C COEFFICIENTS OF ZEROS APE 

AA1=. 101 2235 362 Q03762600 
A A 2=. 2223810344533744700 
AA3=, 3137066458778872900 
AA4=. 362683783378361 98D0 

C XG CONTAINS ALL THE 8 ZEROS AND WT CONTAINS ALL THE 
C COEFFICIENTS. ALL APE SHIFTED TO RANGE 0 TO 1. 

WT(1 ) = AAl=fc0.5D0 
WT{ 2)=AA2*0. 500 
WT(3)=AA3*n.5D0 
WT(4)=AA4*0,500 
WT( 5)=WT( 4) 

WT(6 )=WT( 3) 

WT(7)=WT(2) 

WT( 8)=WT( 1 ) 

A=0.500 

XG( 1 ) =A + XX1’)'C. 5D0 
XG( 2) = A+XX2*C. 5D0 
, XG(3)=A+XX3+0.5D0 
XG(4) = A + XX4>»C.500 
XGI 5 )=A-XX4*0. 500 
XG(6) = A-XX3*-0.5D0 
XG(7)=A-XX2-^C. 500 
XG(8)=A-XX1*0.500 
DO 100 1=1,8 
00 100 J=l,8 
100 WTWI I, J)=WT( I)*WT(J) 

DO 200 1=1,8 
X = XG( I ) 

DO 200 J=l,8 
Y=XG( J) 

CALL FORMK(STF,X,Y) 

RALF = 1 .OCO/QQIX ) 

DO 200 K=l,24 
DO 200 M=l,24 

200 STK(K,M) = STK(K,M)+WTW( J, I ) *S TF { K , M ) ♦L^R ALF 
DO 300 1=1,24 
00 300 J=l,24 

STK( I , J) = (STK( I, J)+STK(J, I ) ) *0.500 
3CC STK( J, I )=STK (I , J) 
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IF(KASF.FQ.O) RFTURN 
CALL TFCRMdlNVfTINVT ) 

CALL MUl T( TINVT,STK,STIFF1 ,.?4,24,2^) 

CALL mulT( STIFF l ,TTNV»STIFFF,?4,24,24) 

on 600 1=1,24 

on 60C J=l,24 

ST I FFF( I , J ) = 0. 500 * (STTFFF( I,J) + STIFFF(J,in 
STIFFF(J,I)=STIFFF(I,J> 

600 CONTINUF 

IF(KASE.EC.2 ) RFTURN 

CALL RIGID ( STl FFFrSTFMODI 

RETURN 

END 



SUBROUTINE FOR MK ( S TF , X , Y ) 

C FCRYK FPRYS THE ELEMENT STIFFNESS MATRIX. 
implicit realms (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

DIMENSION STF(2A,24) , WORK 1 (6,24 ) , OFL ( 6 , 6 I , R ( 6 , 24 ) 
CALL DEL1(DEL,X) 

CALL PST(P,X,Y) 

DO 200 1=1,6 
on 200 J=l,?4 
WORKl ( I ,J 1 = 0.000 
00 200 K=l,6 

200 WORKl ( I ,J)=WCPK1( I,J)+DEL(I,K)*P(K,J) 

DO 300 1=1,24 
DO 300 J=l,24 
STF( I , J)=0.000 
on 300 K=l,6 

300 STF( I , J)=STF ( I , J)+P(K, I )*W0RK1 (K, J I 
RETURN 
END 



SUBROUTINE TFOR M ( T I NV , T I NV T ) 

C TFORM FORMS THE TRANSFORMATION MATRIX AND ITS INVERSE. 
IMPLICIT REAL*8 (A-H,0-Z> 

REALMS K1,K2,K3,L 

OIMFNSION TINVT( 24,24) ,TINV( 24, ’4) 

COMMON/Pl /E,H1 ,H2,POIS,RA01,RA02,L , PHI 1 , PH I 2 , ALPHA, 
#THETA,KASE 

C INITI ALIZE T INVT MATP IX 
DO 350 1=1,24 
DO 350 J=l,24 
350 TINVT( I , J)=O.ODO 

IF( KASE.E0.2 ) ALPHA=WIOTH 
IF( KASE.EQ.2 ) RAD1=1.000 
IF(KASE.EQ.2 ) RA02=1.000 
TINVT( 1, 1)=1.000 
TINVT( 1, 2)=-1.0D0 

TINVT ( I , 3)=-] .000 
TINVT( 1, 4)=1.000 
TINVT( 2, 5)=1.0n0 
TINVT( 2, 6)=-1.0D0 
TINVT( 2, 7)=-1.0D0 
TINVT( 2, 8)=1.0D0 
TINVT( 3 , q)=4.000 

tinvt( :»,io)=-6.odo 
TINVT( ?,12)=2.0D0 
TINVT( 3,13)=-6.0D0 
TINVT( 3,14)=q.0D0 
TINVT( 3,16)=-3.000 
TINVT( 3,21)=2.0D0 
TINVT( 3,22)=-3.0D0 
TINVT( 3, 241=1.000 
TINVT( 4 , q)=2.000 *L 
TINVT( 4,10)=“3.0D0 *L 
TINVT( 4,121=1.000 *L 
TINVT( 4 ,1 3) =-6^. ODO *L 
TINVT( 4,14)=6.0D0 *L 
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TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 



4 


,16) 


=-2.000 


4 


,17) 


=2.000 


4 


,18) 


=-3.000 


4 


,20) 


= 1.000 




, 9) 


=2. 000 




,10) 


=-4.000 




,11) 


=2.000 


5 


,18) 


=-3. 000 


5 




= 6.000 


6 


,15) 


=-3.000 


5 


♦ 21) 


=1.000 


■5 


,22 ) 


=-2.000 


c 


,28) 


=1. 000 


h 


♦ 9) 


=1.000 


6 


,10) 


=-2.000 


6 


,11) 


= 1. 000 




,18) 


=-2.000 


6 


,14) 


=4.000 


6 


,15) 


=-2.000 


6 


,17) 


=1.000 


6 


,18) 


=-2.000 


6 


, 19) 


=1. 000 


7 


, 1) 


=-1.000 


7 


, 2) 


=1.000 


R 


, 5) 


=-1.000 


P 


, 6) 


=1.000 


Q 


♦ 9) 


=-4.000 


9 


,10) 


=6. 000 


9 


,12) 


=-2.000 


Q 


♦ 13) 


=6. 000 


Q 


,14) 


=-9.000 


q 


,16) 


=3. 000 


10 


, 9) 


=2. 000 


10 


,10) 


=-3.000 


10 


,12) 


=1.000 


10 


,13) 


=-2.000 


10 


,14) 


=3.000 


10 


,16) 


=-1.000 


11 


, 9) 


=-2.000 


11 


,10) 


=4.000 


1 1 


,11) 


=-2. 000 


11 


,13) 


=3.000 


11 


,14) 


=-6.000 


1 1 


,15) 


=3. ODO 


12 


, 9) 


=1.000 


1? 


,10) 


=-2. ODO 


12 


,11) 


=1. 000 


12 


,18) 


=-1.000 


12 


,14) 


=2. 000 


12 


,15) 


=-l. ODO 


1? 


, 1) 


=1.000 


14 


, 5) 


=1. ODO 


15 


, 9) 


=4. ODO 


15 


,10) 


=-6.000 


15 


,18) 


=-6. ODO 


15 


,14) 


=9.000 


16 


, 9) 


=-2.000 


16 


,10) 


=3. 000 


16 


, 13) 


=2.000 


16 


,14) 


=-3.000 


17 


, 9) 


=-2. ODO 


17 


, 10) 


=2.000 


17 


,13) 


=3.000 


17 


,14) 


=-3. ODO 


IB 


, 9) 


=1.000 


18 


,10) 


=-1.000 


18 


,13) 


=-1.000 


18 


,14) 


= 1.000 


19 


, 1) 


=-1.000 


19 


, 8) 


=1.000 


20 


, 5) 


=-1.000 


20 


, 7) 


=1. ODO 



*L 

*L 

*ALPHA*RAD1 
*ALPHA*RAD1 
*ALPHA*RA01 
*ALPHA*RAD1 
*ALPHA*RAD1 
*ALPHA*RAD1 
*ALPHA’^RA01 
*ALPHA*RAni 
*ALPHA*RAOl 
*L*ALPHA*RAD1 
^L*ALPHA*RAD1 
•^'L'^Al PHA*RA01 
TL*ALPHA-*RAD1 
*L*ALPHA«RAD1 
*L*ALPHA*RAD1 
*ALPHA'^RAD1 
*L*ALPHA*RAD1 
♦LTALPHA-^RAni 



*-L 

*-L 

*L 

*ALPHA*RAO? 

*ALPHA*RA02 

*ALPHA*RAD? 

*ALPHA«RAD2 

*ALPHA’“RAD? 

*ALPHA*RAD? 

*L*ALPHA«P A02 

*L" ALPHA*RAH2 

«L*ALPHA’^RAD2 

TL*ALPHA*RA02 

*L*ALPHA’^'RAD? 

*L*ALOHA*RAD2 



*L 

*L 

#L 

*L 

*ALPHA*R A02 
*ALPHA*RAD2 
*ALPHA*RAD? 
*ALPHA*RAn2 

♦ L*ALPHA*RA02 
*L*ALPHA*RA02 

♦ L*ALPHA<'RA02 
*L*ALPHA'^RA02 
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TINVT(?1, q»=-^.ODO 
TINVK ?1 ,10)=6.0D0 
TINVK ?1 ,n) =6.000 
TINVT (21 ,14»=-o.000 
TINVTi 21 ,21) =-2.000 
TINVT( 21 , 22) =3. ODO 
TINVT(22, Q)=-2.000 
TINVT( 22,10)=3. 000 
TINVT( 22, 13)=4.000 
TINVT (22, 14)=-6. 000 
TINVT(22,17)=-2.000 
TINVT( 22, ia) = 3.000 
TINVT(23, 9)=2.000 
TINVT(23,10)=-2.000 
TINVT(2^,13)=-3.000 
TINVT (23 ,14) =3.000 
TIMVT( 23,21)=1.000 
TINVT( 23,22)=-1.000 
TINVT (24, 9) =1.000 
TINVT(24,10)=-1.000 
TINVT( 24,13)=-2.0D0 
TINVT(24,14)=2.0no 
TINVT( 24, 17)=1.0D0 
TINVT(24,18)=-1.0n0 
DO 375 1=4,22,6 
00 375 J=l,24 
375 TINVT( I , J)=-TINVT( I, 
on 500 1=1,24 
on 500 J=l,24 
500 T INV( I , J)=TINVT( J , I ) 
RFTijrn 
END 



♦L 

’^L 

*L 

*L 

■Ml 

^Al.PHA’^'RADl 
*ALPHA*RA01 
*ALPHA'^R AOl 
*ALPHA^RAD1 
^ALPHA-^RAOi 
*ALPHA*RA01 
*L'*ALPHA*R ADI 
*L* ALPHA^RAOl 
♦L^ALPHA’^RAOl 
*L* ALPHA=RA01 
•^L^ALPHA'^RAOl 
*L*ALPHA*RA01 



J ) 



SURRnUTiNE DEL1(DEL,X) 

C DELI FnRh^S the CONSTITUTIVE LAW MATRIX. 

IMPLICIT REAL’^8 (A-H,n-Z) 

REAL*8 K1,K2,K3,L 
DIMENSION 0EL(6,6) 

CnMMON/P,l/E,Hl,H2,PniS,RADl,RAD2,L,PHI I , PH I 2 , ALPHA , 
MTHETA , KASE 
H( X)=H1+(H2-H1 )-*X 
Ki=(H(x)*F)/( 1 .ono-poi S**2 ) 

K2 = K1* PCI S 

K3=(K1*( l.OnC-POIS) ) /2.0D0 

oi = (H(x )**3*P)/ ( 1, 2oi*( 1 .ono-pn is**?)) 

n->=ni*pri s 

03=(01*(1.0DC-PniS) ) /2.0D0 
on ICO 1=1,6 
on 100 j=i,6 
100 DEL ( I, J)=c.000 
DEL( 1 , 1 )=K1 
OEL( 1 ,2)=K2 
OEL ( 2, 1 ) = K2 
0EL(2 ,2)=K1 
OEL( 3,3)=K3 
0FL(4,4)=01 
0EL(4,5)=02 
nEL( 5,4)=02 
DEL(5,5)=01 
0EL(6,6)=03 
RETURN 
END 



FUNCTION P 1C1(X,Y) 

IMPLICIT REAL*R (A-H,0-Z) 

REALT-8 KI ,K2,K3,L 

COMMON/R! /E,H1 ,H2 , POI S ,RAD1 ,RA02,L , PHI 1 , PH I 2 , AL PHA, 
<(THETA, KASE 
Plot =Y/L 
RETURN 
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entry P 1C2(X,Y> 

pio? =i.oon/L 

RETURN 

ENTRY P lOP(X.Y) 

P109 =Q1(X )*x**3*Y**3 

RETURN 

ENTRY P 110(X,Y) 

PllO =Q1(X )*x**3*Y**2 

RETURN 

ENTRY P 111(X,Y) 

Pill =Q1(X )*X**3*Y 

RETURN 

ENTRY P 112(X,Y) 

P112 =G1(X >*X*’*'3 

RETURN 

ENTRY P 113(X,Y) 

P113 =01 (X )^X**2*Y**3 

RETURN 

ENTRY P IIA(X.Y) 

P114 =01 (X )*x**Z*y**^2 



RETURN 



ENTRY 

P115 


P 115(X 
= 01 (X 


*Y 


RETURN 

ENTRY 


P 116(X 


tY ) 


P116 

RETURN 


= 01 (X 


I*X**2 


ENTRY 


P 117(X 


tY ) 


P117 

RETURN 


= 01 (X 


) ♦X'*Y**3 


ENTRY 


P 1181 X 


tY ) 


P118 

RETURN 


= 01 (X 


)^X*Y**2 


entry 


P 119(X 


tY > 


PllR 

RETURN 


= 01 (X 


)*X*Y 


ENTRY 


P 120IX 


;^x 


P120 

RETURN 


=01 (X 




ENTRY 


P 12KX 


tY ) 


P121 

RETURN 


=cnx 


)*Y**3 


ENTRY 


P 122IX 


tY ) 


P122 

RETURN 


= 01 IX 


)*Y^*2 


ENTRY 


P 123IX 


tY ) 


P123 

RETURN 


= 01 (X 


)*Y 


ENTOY 


P 124IX 


tY) 


P124 

RETURN 

END 


=01 IX 


) 


EUNfTICN 


p 2C1 I X 


tY) 


IMPLICIT 


R EAL*8 


I A-H,0-Z ) 



REALMS K1 ,K2 fK3,L 
COMMON/Bl/EfHl ,H2tPniS,RADl 
#THETA,KASE 



P201 

RETURN 


= ICCSPHI I X 


)*OIX 


entry 


P 202IX,Y) 




P202 

RETURN 


=ICCSPHI IX 


)*OIX 


ENTRY 

P203 


P 203IX,Y) 
=ICCSPHIIX 


)*OIX 


return 

ENTRY 


P 204IXtY) 




P204 

return 


= ICCSPHHX 


)*Q{X 


ENTRY 

P205 


P 2051 X.Y ) 
=OOIX )*X 





, RAD2tL,PHIl , PH I 2, ALPHA, 
) )*X^Y 

) )*X 

) ) ^Y 

) ) 
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RETURN 

ENTRY 


P ?07(X,Y) 






P207 

RETURN 


=CQ(X ) 






ENTRY 


P 2C9(X,Y) 




) x#*3^Y** 3 


P20<5 

RETURN 


=(SINPHI(X 


)*Q(X 




ENTRY 


P 210( X,Y) 






P210 

RETURN 


=( S INPHI (X 


) >!<0(X 


) ) *X **3*Y*” 2 


ENTRY 


P 21K X,Y) 






P211 

RETURN 


=(SINPHI(X 


) *Q(X 


) ) AX’!f*3*Y 


ENTRY 


P 212( X,Y ) 






P212 

RETURN 


=( S INPHI ( X 


)*Q(X 


) )>iX'^*3 


ENTRY 


P 213(X,Y) 






P2] 3 
RETURN 


=(SINPHI(X 




) ) *X'^*2^Y**3 


ENTRY 


P 214(X,Y) 






P23 4 
RETURN 


= (S INPHKX 


) *0(X 


) )*X-‘*?<'Y**2 


ENTRY 

'»215 

RETURN 


P 215(X,Y) 
= ( S INPHI (X 


)TQ(X 


) ) *X**2’^Y 


ENTRY 


P 216(X,YI 






P216 

RETURN 


=(SINPHI (X 


)*Q(X 


) 


ENTRY 


P 217(X,Y) 






P217 

RETURN 


= ( S INPHI (X 


) *0(X 


) Jx-X*YT*3 


ENTRY 


P 218( X,Y ) 






P21 0 

RETURN 


=(S INPHI (X 


)*Q(X 


) )^x*Y*t2 


ENTRY 


P :’19(X,Y) 






o 

RETURN 


=(SINPHI(X 


)*0(X 


) )^X*Y 


ENTRY 


P 220(X,Y) 






P220 

RETURN 


= ( S INPHI ( X 


)’*Q(X 


) )‘X 


ENTRY 


P 22K X,Y) 






P221 

RETURN 


= (SINPHI (X 


)^Q(X 


) )*y**3 


entry 


P 2?2(X,Y) 






P?2R 

RETURN 


=(SINPHI(X 


)*0(X 


) ) ’^Y**? 


ENTRY 


P 223(X,Y) 






P223 

RETURN 


^(SINPHKX 


) +Q(X 


) )*Y 


ENTRY 


P 224(X,Y) 






RETURN 

END 


=(SINPHI (X 


)*Q(X 


) ) 



FUNCTICN P 3C1(X,Y) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*B K1,K2,K3,L 

CQMM0N/B1/E,H1 ,H2»PniS»RADl,RAD2,L, PHIl ,PHI2» ALPHA, 
«THETA, KASF 



P301 

RETURN 


=CQ(X )+X 






ENTRY 


P 303( X , Y 1 






P303 

RETURN 


=CQ(X ) 






ENTRY 


P 305(X,Y) 






P305 

RETURN 


= V/L-(0(X 


)/L)*DR(X 




entry 


P 306(X,Y) 






P306 

RETURN 


=1.000/L-(0(X 


)/L)*nR(X 


> *X 


ENTRY 


P 307( X, Y) 
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P^07 =(-Q(X )/L)*DR(X )*Y 

RETURN 

ENTRY P 30R(X,Y) 

P308 =(-C(X l/L)*OR(X ) 

RETURN 

END 



FUNCTION P 4C1 ( X,Y» 
implicit REAL*8 (A-H,0-Z) 

REAL*8 KI,K?,K3,L 

COMMON /Rl/E, HI , H2 , PO I S t R AD 1 , RAD? , L , ph I 1 ,PHI2 
«THETA , KASE 

PAOl =-(Gl(X )**2/L )*DRl( X I'^'X^Y + IQK X 

RETURN 

entry P 402(X,Y» 

PA02 =-(Cl(X )**2/L)*OPl (X )*X +(Ol(X 

RETURN 

entry P A03(X,Y) 

P403 =-(Cl(X )-^*2/L)*DRl (X >^Y 

RETURN 

entry p A0A(X,Y) 

P40A =-(Ql(X )*^T2/L)*nRl (X > 

RETURN 

entry p 40P( X, Y ) 

P409 =(-6,ODO/L**2 ) *X=»Y*"'3 

RETURN 

entry p AlO(X,Y) 

PAID =(-6,0n0/L*’^2 )"'X*Y*^2 

RETURN 

entry P A11(X,Y) 

PAll =(-6.0D0/L**2 )*X*Y 

RETURN 

ENTRY P A12(X,Y> 

PA12 =(-6.0D0/L**2 )*X 

RETURN 

entry P A13(X,Y> 

PA13 =( -2.nnO/L**2 )*Y^»3 

RETURN 

ENTRY P 41A(X,Y) 

PAIA =(-2.0n0/L**2 )<'Y**2 

RETURN 

ENTRY P A15( X ,Y » 

PA15 =(-2.nnO/L*’>‘2 ) *Y 

RETURN 

entry P 416( X,Y ) 

PA16 = (-2.OD0/LT + 2 ) 

RETURN 

ENTRY P A90(X,Y) 

P4qO =COSPHI(X »*Q(X )/L 

RETURN 

entry P A91(X,Y) 

R491 =CO(X )**?. 

RETURN 

END 



FUNCTION P 5C1 ( X, Y) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 K1,K2,K3,L 

COMMON/Bl/E.Hl ,H2,POIS,RADl,RAn?,L,PHIl ,PHI2 
«THETA,KASE 



P501 


=COSPHI (X 


)+0(X 


» -^01 ( X 


) *X'-^ Y 


RETURN 

entry 

P5D2 


P ■502(X,Y) 
=COSPHI (X 


) + 0(X 


)^01 (X 


) «X 


RETURN 

ENTRY 

P503 


P 503(X,Y) 
=CnSPHI (X 


)*Q(X 


)TQ1 (X 


) *Y 


RETURN 

ENTRY 

P50A 


P 50A(X,Y) 
=COSPHI (X 


)*Q(X 


)*Q1 (X 


» 



ALPHA, 

)/L)*Y 

)/L) 



ALPHA, 
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RETURN 

ENTRY 

P505 

RETURN 


P 505(X,Y) 

=Q(X »+QQ(X )*X 


ENTRY 

P507 

RETURN 


P 507(X,Y) 

=G(X )*00(X ) 


entry 

P509 

i^X**T*Y 


P 509(X,Y) 

= -3.0DO*P490(X,Y > •<'X**2'*'Y**3-6. 0D0*P491 ( X, Y) + 


RETURN 

ENTRY 

P510 

MV 

RETURN 


P 510(X,Y) 

= -3,0D0*P490 ( X , Y » ^ X** 2* Y** 2 -2 . 0D0*P491 ( X ,Y ) * 


entry 

0511 

RETURN 


P 51K X,Y> 

= -3.0nO*P490( X, YJ «X**o «Y 


ENTRY 

P512 

RETURN 


P 512( X,Y) 

=-3.0D0*P490(X, Y»*X**2 


ENTRY 

P513 

H*y 


P 513( X,Y) 

= -2.0D0*P49O( X, Y) XX* Y’'- *'3-6. 000*P4'31 ( X , Y ) * X* * 2 


RETURN 

ENTRY 

P51A 

RETURN 


P 514( X , Y ) 

= -2.0D0*P490(X,Y) 'X *Y* *2-2 . 000* P49 1 (X,Y )*X**2 


ENTRY 

P515 

RETURN 


P 515 (X,Y ) 

=-2.ODO*P4g0( X,Y) *X*Y 


entry 

P516 

RETURN 


P 516 ( X, Y » 
=-0,0D0*P490( X,Y) '^X 


ENTRY 

P517 

RETURN 


P 517(X,Y) 

= -P4Q0( X, Y)'*Y **3-6,000 XP4Q1 {X,Y»*X*Y 


entry 

P518 

RETURN 


P 518 ( X,Y ) 

=-P490( X, Y»*Y**2-2. 000*P4Ol (X , Y )*X 


entry 

P519 

RETURN 


P 519(X,Y) 
=-P490( X, Y) " Y 


ENTRY 

P520 

RETURN 


P 520(X,Y) 
=-D490( X, Y) 


ENTRY 

P521 

RETURN 


P 521 ( X,Y ) 
=-6,000*P491 ( X, Y) «Y 


entry 

P522 

RETURN 


P 522 (X,Y» 
=-2.0D0*P491 ( X,Y) 


pNTRY 

P5R0 

return 


P 590( X,Y) 

= (Q(X )**2/U*0R(X ) 


entry 

P5R1 

RETURN 


P 591 ( X,Y ) 

= (0(X )*QQ(X )/U*nR(X ) 


ENTRY 

P592 

RETURN 

END 


P 592 ( X, Y ) 
=QQ(X )/L 


FUNCT ION 


P 6C1 ( X, Y) 



IMPLICIT REAL*8 (A-H,0-Z» 

RFAL<-8 K1,K2,K3,L 

C0MM0N/B1/E,H1,H2, POI S,RAD1 ,RAD2,L . PHI 1 ,PHI2 , ALPHAt 
«THETA,KASF 



P601 

RETURN 


=(GC(X )*01(X ))*X 


ENTRY 


P 603(X,Y) 
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P603 =CQ<X )’*'Q1(X ) 

RETURN 

ENTRY P 605(X,Y) 

P60*^ =-2.0D0*P590( X,Y) *X^Y+ ( Q< X )/L)*Y 

RETURN 

ENTRY P 606(X,Y) 

P606 =-?,0n0*P590(X,Y) *X +(0(X )/U 

RETURN 

ENTRY P 607(X,Y) 

P607 = -2. 0D0*P590( X,Y) *Y 

RETURN 

•^NTPY P 60R ( X,Y » 

P608 = -2.000*P590(X,Y) 

RET(JRN 

ENTRY P 609(X,Y) 

P60R =6, 0D0*P5P1 ( X ♦ Y»*X«*3<'Y+>»=2-18.0D0*P592 ( X, Y ) * 

#X'**2*Y**2 
RETURN 

ENTRY P 810<X,Y» 

P610 =A.000*P591(X,Y)*X**3*Y - 1 2. ODO^ P 592 ( X , Y ) * 

#X**2*Y 
RETURN 

entry P 51i(X,Y> 

P611 =2.000*P5 9UX,Y)TX'**3-6.0DO*P592< X, Y»*X**2 

RETURN 

entry P 613( X,Y ) 

P6J.3 =5, 0D0*P5 91 (X,Y ) *X«t 2>:- YT Jt2- 1 2 . 000 *P 592 ( X, Y ) ^ 

«X*Y**2 
RETURN 

ENTRY P 614( X,Y ) 

P61A =4. CD0*P5 91 ( X , Y ) * X*«2^ Y-8. 0 00+ P592 ( X , Y ) -^X^Y 

RETURN 

ENTRY P 615( X,Y) 

P615 =2. 000*P591 ( X,Y>*X*-^2 - 4. 000* P59? ( X , Y ) *X 

RETURN 

ENTRY P 517(X,Y) 

P617 =6. CD0*P5Q1 ( X, YJ + X^Y*- 2-6. 000* P5P2 ( X » Y ) *Y**2 

RETURN 

ENTRY P 618( X,Y» 

P618 =4. 0D0*P591 ( X, Y)*X *Y - 4 . 000* P59? ( X , Y ) *Y 

RETURN 

ENTRY P 619< X,Y) 

P619 =2. CDO*P591 ( X, Y)*X - 2 . 000* P5 92 ( X ♦ Y t 

RETURN 

ENTRY P 621 < X,Y) 

P621 =6. C00*P591 ( X, Y»*Y**2 

RETURN 

ENTRY P 622 ( X, Y) 

P622 =4. 0D0*P59H X,Y)*Y 

RETURN 

ENTRY P 6?3( X,Y) 

P623 =2.000*P59I (X,Y) 

RETURN 

END 



SUBROUTINE PST(P,X,Y) 

C PST CONVERTS THE P'S EROM EUNCTIONS OF X AND Y TO VALUES 

C FOR ELEMENTS QE THE P MATRIX. 

IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 K1,K2,K3,L 
DIMENSION P(6,24) 

P(l, 1)=P101(X,Y) 

P(l, 2)=P1C2(X,Y) 

P(l, 3)=0.0DC 
P(l, 41=0.000 
P(l, 5)=C.CDC 
P(l, 6) = O.OOC 
P(l, 7)=0.000 
P(l, R»=O.ODC 
P(l, 9)=P109(X,Y) 

P( 1 ,10) = P110(X,Y) 
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RETURN 

ENTRY 

P505 

RETURN 


P 505(X,Y) 

=G(X )TQQ(X )*X 


ENTRY 

P507 

RETURN 


P 507(X,Y) 

=G(X I*00(X ) 


entry 

P509 

(^X**^*Y 


P 509(X,Y) 

= -3.0D0*P490 (X , Y 1 ■*X**2'*'Y**3-6. 000*PA91 ( X, YH' 


RETURN 

ENTRY 

P510 

uv **-1 

RETURN 


P 510(X,Y) 

= -3,00 0*P490( X, Y) '«‘X**2*Y**2-2. 0D0*P491 ( X ,Y ) # 


entry 

P511 

RETURN 


P 511( X,Y) 

= -3.0n0*P49 0(X,Y)i'X**'>ftY 


ENTRY 

P512 

RETURN 


P 512( X,Y> 

=-3.0D0*P490( X, Y) *X**2 


ENTRY 

P513 

«*Y 


P 513(X,Y) 

= -2.0n0*P490( X, Y) xX^Y>' *-?-6. 0D0*P4'3l ( X, Y)*X**2 


RETURN 

PNTRY 

P51^ 

RETURN 


P 514( X, Y I 

= -2.0D0*P490(X,Y) 'X *Y * *2-2 . 000* P49 1 (X,Y »*X**2 


entry 

P515 

RETURN 


P 515 < X,Y > 

=-2.0n0*P490( X,Y) *X*Y 


ENTRY 

P516 

RETURN 


P 516 < X,Y » 
=-7.000*P490( X,Y) «X 


ENTRY 

P517 

RETURN 


P 517(X,Y) 

=-P4P0(X,Y)*Y**3-6.0D0*P4Ql(X,Y»*X*Y 


entry 

P518 

RETURN 


P 518(X,Y ) 

=-P490(X, Y)*Y**2-2.000*P4P1 (X,Y)*X 


entry 

P519 

RETURN 


P 519(X,Y» 

= -P490 ( X, Y) Y 


ENTRY 

P520 

RETURN 


P 520(X,Y> 
=-P490( X, Y) 


ENTRY 

P521 

RETURN 


P 521 ( X,Y ) 
=-6.0D0*P491 ( X,Y) »Y 


ENTRY 

P522 

RETURN 


P 522(X,Y) 
=-2.000*P491 ( X,Y> 


ENTRY 

P590 

RETURN 


P 590(X,Y» 

= (Q(X )**2/U*nP(X ) 


ENTRY 

P591 

RETURN 


P 591 (X,Y ) 

= (0(X )*QQ(X )/U*DR(X ) 


entry 

P592 

RETURN 

END 


P 592 ( X,Y ) 
=QO(X )/L 


FUNCTION 

IMPLICIT 


P 6CK X, Y) 

REAL*« (A-H,0-Z) 



RFAL=<'‘^ K1,K2,K3,L 

C0MM0N/Bl/E,Hl»H2t P0IS,RAD1 ,RAD2»L, PHIl ,PHI2 » ALPHA, 
#THETA, KASF 



P601 

RETURN 


=(QQ(X )*01(X ))*X 


ENTRY 


P 603(X,Y) 



69 



P603 

RPTUPN 

ENTRY 

960*^ 

RETURN 

entry 

P606 

RETURN 

ENTRY 

P607 

RETURN 

«=NTPY 

9608 

RETURN 

ENTRY 

P60R 

#X'**2*Y**2 

RETURN 

ENTRY 

P610 

#X**2’^Y 

RETURN 

ENTRY 

P611 

RETURN 

entry 

P613 

«X*Y**2 

RETURN 

entry 

P6I4 

RETURN 

ENTRY 

P615 

RETURN 

ENTRY 

P617 

RETURN 

ENTRY 

P618 

RETURN 

entry 

P619 

RETURN 

ENTRY 

P621 

RETURN 

ENTRY 

P622 

RETURN 

ENTRY 

P623 

RETURN 

END 



=C0(X »*01(X I 

P 605(X,YI 

=-2.OD0*P590( X,Y) ♦X^'Yf (Q(X »/L)*Y 
P 6 06 ( X , Y » 

= -2,000<'P5O0( X,YJ *X +(0(X l/L) 

P 607( X»Y» 

= -2.0DO»P590(X,YI*Y 

P 6 08 ( XtY » 

= -2. 000*P590( X,Y I 

P 609( X, Y » 

=6,000*P5Pl(X,Y»*X**3*Y**2-18.0n0^P592(X,Y)* 



P 610( X, Y» 

=4.000*9591 (X»YI*X**3*Y - 1 2. 0D0*P 592 ( X , Y » * 



P 6 1 i ( X , Y » 

=2.0D0*P591 (X, Y )*X **3-6.000*9592 ( X, Y)*X**2 
P 613( X»Y » 

= 6.000*95 91 (X,Y J*X = *2*Y*)«2-12.n00*P592( X, Y»* 



P 614(X,Y» 

=4. COO* P5 91 (X,Y»*X**2*Y-8.000*P592( X,Y )*X*Y 
P 615( X,Y» 

=2. ODO*P591 ( X,Y)*X**2 - 4. 000* P592 ( X , Y ) *X 

D 617( XtY) 

=6. CD0*P5P1 ( X» Y)*X*Y*" 2-6. 000* P592 ( X , Y ) *Y**2 
P 618( XtY) 

=4. 0D0*P591 ( Xt Y)*X *Y - 4. 000*959? ( X t Y ) *Y 
P 619( X,Y) 

=2. CD0*P591 ( Xt Y)*X - 2. ODO*P5R2 ( X tY) 

P 621 ( XtY) 

= 6. COO* P 5 91 ( X t V)* Y**2 

P 622 (XtY) 

=4. 0D0*P591 (Xt Y)*Y 

P 6?3( X,Y) 

= 2. 0D0*P591 (XtY) 



SUBROUTINE PST(PtXtY) 

C PST CONVERTS THE P'S EROM EUNCTIONS OE X AND Y TO VALUES 

C EQR ELEMENTS OF THE P MATRIX. 

IMPLICIT REAL*e (A-HtO-Z) 

REAL*8 KltK2tK3tL 
DIMENSION P( 6t 24) 

P(lt l)=P101(XtY) 

P(lt 2)=P102(XtY) 

P(lt 3)=0.0DC 
P(lt 4)=0.0D0 
P(lt 5)=C.CDG 
P(lt 6) = C.OOC 
P(lt 7)=0.0D0 
P(lt 8) = O.OOC 
Pdt 9)=P109(XtY) 

P(ltl0)=P110(XtY) 



70 



p( i,in = piii (x,Y) 

P( 1,12)=PH2(X,Y) 
P(ltl3) = Pin(X,Y) 
P( 1 ,14)=P114{X,Y) 
P{1,15»=P115(X,Y) 
P(l »16 ) = Plli'. (X,Y) 
P( lt!7)=PU7(X,Y» 
P(1,18)=P113(X,Y) 
P(1 ,?P)=P11P(X,Y) 
P( \ ,20)=P120(X,Y) 
P(1,21)=P121(X,Y) 
P(1,22)=P122(X,Y» 
P( 1,23)=P123(X,Y) 
P(1,2^)=P12A(X,Y» 
P(2, 1 )=P?01 (X,Y) 
P(2, 2)=P2C2(X,Y) 

P(2, 3»=P203(X,Y» 
P(2, A)=P2C4(X,Y) 
P(2, 5)=P205(X,Y) 

p(2» 6)=c.ono 

P(2, 7)=P2C7(X,Y) 

p( 2 , a)=c.onn 

P(2, P)=P20c>(X, Y > 
P( 2f 10)=P210(X, Y» 
P( 2, 1 H = P211 ( X, Y ) 
P(2,12)=P212(X,Y) 
P(2,13)=P213(X,Y) 
P( 2, 14 ) = P214( X, Y ) 
P(2 ,15 ) = P215 (X,Y ) 
P( 2, 16)=P216(X, Y) 
P(2,17)=P217(X,Y> 
P(2tl P)=P218(X,Y) 
P(2,1P)=P21'3(X,Y) 
P{? ,20 ) = P220(X, Y ) 
P(2 ,21 ) =P221 (X,Y ) 
P(2,22)=P222(X,Y) 
P(2 ,23 )=P223(X,Y ) 
P(2,241=P224(X,Y) 
P(3, n = P301(X,Y) 
p(3, 2)=c.ono 
P(3, 3)=P3C3(X,Y1 
P(3, 4)=C.0D0 
P(3, 5)=P305(X,YI 
P(3, 6»=P3C6(X,Y) 
P(3, 7I=P3C7(X,Y» 
P(^, 8)=P30S(X,Y) 

P(3, P)=0.0DC 
P(3,10)=0.CDC 
P(3 ,11 )=0.000 
P( 3,12»=c.cnc 
P(3,13 ) = n,0D0 
P(3 ,14)=0.0n0 
P(3,15»=C.COO 
P(3 ,16 1=0.000 
P(3,17)=0.0D0 
P( 3,18)=C.C0C 
P( 3 ,10 )=O.ODO 
P(3 ,20) =0.000 
P( 3,21 )=C.C0C 
P(?,22)=0.000 
P( 3,23) =C. COO 
P( 3,24)=C. COC 
P(4, 1)=P401(X,Y) 

P(4, 2)=P402(X,Y) 
P(4, 3)=P403(X,Y) 
P(A, 4)=P404(X,Y) 
P(4, 5)=o.ono 
P(4, 6)=0.000 
P(4, 7)=0.00n 

p<4, 8)=c.ono 

P(4, 9)=P400(X,Y) 
P(4,10)=P410(X,Y) 
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P(4,11)=P411(X,Y) 
P(4,1?)=P412(X,Y) 
P(4,13)=P413(X,Y) 
P(4,14)=P414(X,Y) 
P(4,15 )=P41‘>(X,Y ) 
P(4,16)=P41<S(X,Y) 
P(4, 17)=0.0DC 
P(4,18)=0.0D0 
P(4,19)=0.0DC 
P(4,20)=C.0DC 
P(4,21 )=0.0D0 
P(4,22)=C.CD0 
P(4, 23)=C.00C 
P(4,24)=C.000 
P(5, 1)=P5C1(X,Y) 
P(S, 2)=P502(X,Y) 
PC'), 3)=P503 (X,Y) 
P<5, 4)=P504(X,Y) 
P(5, 5)=P505(X,Y) 
p(5, 6)=o.onn 
P<5, 7)=P5C7(X,Y) 
P(5, 8)=0.0D0 
P(5, 9)=P509(X,Y) 
P< 5,lC)=P5in(X,Y) 
P(5,ll ) = P511 (X,Y» 
f>(5,12)=P512(X,Y) 
P< 5,!3)=P513(X, Y) 
R(5,14)=P51^(X,Y) 
P< ‘5,15)=P515(X,Y) 
PC5, 16) = P5l'S(X, Y) 
P< 5 , 17 ) = P517 (X, Y) 
Pi 5,18)=P51R(X,Y) 
P(5,19)=P519(X,Y) 
P(5 ,?0)=P520(X,Y) 
P( 5,2U=P521 (X,Y) 
P(5,22)=P522(X,Y) 
f’(5»23)=0.0D0 
Pi 5,24)=0.0DC 
P(6, 1)=P601(X,Y) 
P(6, 2)=O.ODO 
P(6, ?)=P6C3(X,Y) 

P(6, 4)=0.000 
P(6, 5)=P6C5(X,Y) 
P(6, 6)=P6C6(X,Y) 

P<6, 7)=P607(X,Y) 
P<6, 8)=P60fl(X,Y) 
P(6, 9)=P609(X,Y) 

P(6,10)=P610(X, Y) 
P<6,11)=P611 (X,Y) 
P(6, 12 ) = 0.0D0 
P(6,13)=P613(X,Y) 
P(6,14)=P614(X,Y) 
P ( , 1 5 ) = P 6 1 5 ( X , Y ) 

P(ft tl6 )=0.0D0 
P(6,17)=P617( X, Y) 
P((S,18) = P618(X,Y) 
P(<S,19)=P(S19(X,Y) 
P( 6,2n)=C.0D0 
P(6,21 ) = P62UX,Y ) 
P(6,22)=P622(X,Y) 
P(6,23)=P623( X,Y) 
P(6,24 )=0.0D0 
RETURN 
END 



SUBROUTINE RIGID ( ST I FFF,STFMOD ) 

C RIGID FORMS THE RIGID BODY AND THE K22ITKT) MATRIX. 
IMPLICIT REAL+8 (A-H,0-7.) 

REAL*8 Kl,K2tK3,L 

DIMENSION TTT( 24,6) tSTFRBI 24t24) ,AKTI24,6) ,TKTINV(6 ,6 ) 
RBI 24, 6) ,ST IFFFI 2 4, 24) , ST F MOD I 24, 24) , TKT ( 6 , 6 ) , A K ( 36 ) , 
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<<AKTT{6,24) 

COMMON /P 1 ,H2 ,Pni S ,RA01 , RAO’ ,L t PHI 1 , PHI2 , ALPHA, 

«THFTA, KASF 

COWMON/P?/CPHI 1 ,CPHI2,C0,CI,C2 
FQUIVALFNCE ( TKT( 1 , 1 ) , AK( 1 ) I 
F7. 1 X I =L'-rSGRT ( 1 . ODO- < CO + C1 ^X+C 2 2 ) * 2 I 

CO=CPHI 1 

IPIKAS F.GF. 1 ) 01=0.000 
IFIKASP.GF.l I 0,2 = 0.000 
IF (KASF.GF. 1 ) CPHI1=0.000 
IFIKASF.EO.^ I CPHI l=(RA02-RA01 ) /L 
IFIKASF.GF.l I CPHI2=CPHIl 
SPHI 1=0 SORT! 1.0n0-CPHIl**2) 

SPHI2=0S0RT( l.OOO-CPHI 2**2 ) 

C NIJMFRICALLY F'/ALUATF FUNCTION Z. 

A = 0. 500 

C = ,4R01 449 24P76R 1200 

7 = . 5061426 51 451 881 200-1 IF Z< A + 0 » ^FZ I A-C ) ) 

C = . 29 R2’2”5706 8133700 

Z=7 + . nil 90 =172? 66 872400 *<FZ<A + 0)+FZ(A-C)) 

C=. 2627662043581644900 

Z=7+. 1 56O5232292894264n0*(FZ( A+C)^F7(A-C) ) 

C=. 91 7 17321 247824900-1 

Z=Z+.l®1341«916aai 800900* ( F7( A+C )+FZ ( A-C ) ) 

TFIK ASF.FQ. 1 I Z = L 
IF<KA5F.E0.2 ) 7=L 
IFIKASF.FO.R ) Z=L*SPHI1 

0 INITIALIZE RIGIC fiOOY MATRIX ANO K2?(TKT) MATRIX. 

no 1 1 = 1 ,6 
on 1 J=l,6 



TKT ( I , J » = 


O.ODO 


00 2 1=1, 


24 


00 2 J=l, 


6 


RBI I , J )=0 


.000 


RBI 1,11= 


CPHI 1 


R8I 1,3)= 


-SPHI 1 


RBI 1,5)= 


Z*RBI 1, 1 l-RAOnRBI 1, 3) 


RBI ?,2)= 


1.000 


RBI 2,4)= 


-Z 


RBI 2,6)= 


R AOl 


RBI 2,1)= 


SPHIl 


RBI 2,3)= 


RBI 1, 1 ) 


RBI 3,5)= 


ZYRBI3, 1 )-P.A01*PBI3, 2) 


RBI 4,1)= 


-II RBI 1 ,1 ) ZRBI3 ,1 ) )*C1 ) /L 


RBI 4,^)= 


Cl/L 


RBI 4,5)= 


-RBI3, 1 )*'-2 + Z^^RBI4, 1 ) -CO* RBI 1 , 1 ) - I R AOl *C 1 ) /L 


RBI 5,2)= 


RBI 3, 1) ZRAOl 


RBI 5,41= 


-RBI 3, 5) /RAOl 


RBI 6,2)= 


RBI4, 1 ) /RADI 


RBI 6,4)= 


-RBI 4, 5) /RAOl 


RBI 7, 1 ) = 


CPHI2 


RBI 7,2)= 


-SPHI2 


RBI 7,5)= 


-RBI 7, 3) *RA02 


RBI R,2)= 


1 .000 


RBI 8,6)= 


RA02 


RBI 0,1) = 


SPHI2 


RBI 0,3)= 


CPHI2 


RBI 9,*') = 


-RA02*RBI9,3) 


RBI 10, 1 ) = 


-I IRBI 7, 1 ) /RBI 9, 1 ) )*ICl+2. 000*C2) ) /L 


RBI 10,3) = 


I Cli-2.0D0*C2)/L 


RBI 10,5) = 


-RBIO,1)**2-ICO+C1+C2)*RBI9,3)-IRA02*IC1+2.0 


RBI 1 1 ,2 ) = 


RBI9, 1 ) /RA02 


RBI 11 ,4) = 


RBI9, 3) 


RBI 12,2) = 


RBI 10,1) /RAD2 


RBI 12,4) = 


-RRI10,5)/PA02 


RBI 13, 1 ) = 


RBI 7, 1 ) ♦OCOSI ALPHA) 


RBI 13,2) = 


RBI 7, 1 ) ♦nSINI ALPHA) 


RBI13,2)= 


RBI7,3) 


RBI 13,4) = 


RA02. + DS INI ALPHA) *RBI 2,2 ) 


RBI 13, 5) = 


R A02+RBI 9, 1 )*OCOSI ALPHA) 


RBI 14, 1 ) = 


-OSINI ALPHA) 


RBI 14, 2) = 


OCOSI ALPHA ) 
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RAD2 

RR(9,1)*RB(14, 
RR( 15»2)= -RR( 9» 1)*RB{ 
Rfl(l'^,'?)= RR(7,1) 

RR(15,^)= RR( 13,2 )*RA02 
RP(15,5)= -RR( 13,1)*’RAD2 
RR(16,1)= RB( 10, 1 J^RRC M 
RR(1^,2)= -RB(10,])*RB{14 
RB(10,3) 
RR(14,1 >»RR(10 
RR( 14,2 )’^RR{ 10 
RR( 9, 1 )*RB( 14, 
RB( 15,1 J/RA02 
RB( 13,1 ) 

RD( 13,2 ) 
RB(10,1 »t'RR(14 
RB( 16, 1 » /RA02 
RR( 14, 2 )*RB ( 12 
RB(1R,5)= -RR( 14,1 )*RB(12 
RB(19,1)= PB( 1, 1)*RR( 14, 
RB(19,2)= -RB( 3, 3)-«^RB{ 14, 

‘ RP(1,3) 

RR( 14,1 )*PB( 1 , 
RB( 14,2 )*RR{ 1, 
RR{ 14,1 ) 

RB( 14,2 ) 

RR(2 ,4) *RB(14, 
RR(20,'=)= -RB( 2,4)+RB( 14, 
RB(?0,6)= RR(2,6) 

RB(21,1)= RB(?,1)*RB(14, 

RB{21,?)= -RB( 3,1)*RB( 14, 
RB(21,3)= RR(1,1) 

R8U4,1 )*PB (3, 
RR( 14, 2)*RB( 3, 
RB(4, 1 ) *RB( 14, 
RB(22,2)= -RB(4, 1 ) >!'RB{ 14, 
RR(22,'’)= RB(4,3) 

RP( 14,])’>-RB{4, 
RB(14,2)*RB(4, 
R8( 5, 2) *PP( 14, 
RB(5,2 )*RR{ 14, 
RB{ 14,2)’^RR( 5, 
RR(23,5)= -PR( 14, n*PR { 5, 
RR(24,1)= RB(4,1)*RB{14, 

RB(24,2)= RB(22,1)/RA01 

R3(24,4)= RP( 14,2)*PB(6, 
RB(24,5)= -RR( 14»U*RB{6, 
DO 3 1=4,22,6 
no 3 J=l,6 

3 RB(I , J)=-RB( I,J) 
1F(KASF.E0.0 ) GO TO 5 
IP(KASE.GT.l ) NDIM=5 
DO 4 1=1,24 

DO 4 J=l,3 

JL=?+J 

JR=3+J 

4 RB( I, JL ) = RB( I, JR ) 
IF(KASE.EQ.1 ) ND1M=4 
GO TO 7 

5 N0IM=6 

7 CONTINUE 
DO 10 1=1,24 
DO 10 J=1,NDIM 
AKT{ I , J)=0.000 
no 10 K=l,24 

AKT( I, J)=AKT(I,J)+STIFFF( 
10 CONTINUE 

DO 20 I=1,NDIM 
DO 20 J=1,NDIM 
TKT( I , J)=C.000 
DO 20 K=l,24 

RR(K,I )=RB( I ,K) TRANSPOSE 



RBI 14, 6 ) = 
RBI15,n = 



R8( 16, 3)= 
RR(16,4)= 
RRI 16, 5)= 
RRI 17,1 ) = 
RRI 17,2 )= 
RRI 17,4)= 
RBI17,5)= 
RRI18,1)= 
RRI 1R,2»= 
RBI1R,4)= 



RBI 10, 3) = 
RBI 1°, 4) = 
RBI19,5»= 
RRI 20, 1 ) = 
RBI20, 2 ) = 
RBI20,4)= 



RRI 21 ,4) = 
RRI 21,5)= 
RBI22,1 )= 



RRI2?,4)= 
RBI 22, 5) = 
R8I23, 1 ) = 
RBI23, 2 ) = 
RRI23,4)= 



2 ) 
1 ) 



, 2 ) 
tl ) 



,5) 
f 5) 

1 ) /RAD2 



, 1 ) /RAD? 

t4) 

,4) 

2 ) 

1 ) 

5) 

5 ) 



2 ) 

1 ) 

2 ) 

1 ) 

5 ) 

5 ) 

2 ) 

1 ) 

5 ) 

5 ) 

1 ) 

2 ) 

4 ) 

4 ) 

1 )/RA01 

4 ) 

4 ) 



I ,K )«RRIK, J ) 
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TKT( I , J)=TKT ( I , J)+R8(K» I )*AKT(K , J) 

20 CGMTINUF 

IF(NDI^'.EC.6 ) GO TO ?5 
on ?1 1=1,5 
I<;T=5*I+1 
no 21 J=IST,35 
JP=J+1 

21 AK(J»=AK(JP) 

IF(NniM.EQ,5 ) GO TO 25 
no 22 1=1,5 
IST=4*1+1 

no 22 J=IST,35 
JP=J+l 

22 AK(J)=AK(JP) 

25 CONTINUF 

CALL INVERT ( NO I M , 1, 00-1 2 , TK T , T K T I NV , KER , NO I M ) 
IF(KFR.EQ,2) WRITF (6, 1000) 

CALL MULT (4KT,TKTINV,TTT,2A,NDIM,N0IM) 

CALL MATT! AKT, AKTT,?4,N0IM) 

CALL MUl T (TTT, AKTT,STFR8,2A,NDIM,24) 

DO 50 1=1,24 
on 50 J=l,24 
STFMODI I , J)=0.000 

STFMOOn ,J)=STIFFF{I,J)-<;TcrB(I,J) 

50 CONTINUE 

Or> 60 1 = 1,24 
on 6C J=l,24 

STFMOOI I , J ) = (STFMOD( I , J ) +S TF MOD ( J , I ) H'O. 500 
STFMPOIJ, I) = STFM00n,J) 

60 CONTINUE 

1000 FORMAT (/,?3H MATRIX K22 IS SINGULAR) 

RETURN 

END 



SUBROUTINE M ATT ( A, AT ,N1 , N2 ) 

C MATT IS AN ONLINE MATRIX TRANSPOSE ROUTINE. 
IMPLICIT REAL*8( A-H,n-Z ) 

DIMENSION A(N1,N2),AT(N2,N1) 

00 ICO 1=1, N2 
00 100 J=1,N1 
100 AT{ I , J) =A( J, I ) 

RETURN 

END 



SUBROUTINE I NV ERT ( N , E P , A , X , KER , NACT ) 
IMPLICIT REAL*8( A-H,0-Z ) 

DIMENSION A (NACT, NACT) ,X (NACT, NACT ) 

00 1 I=1,N 
00 1 J=1,N 

1 X(I,J)=O.ODO 
DO 2 K=1,N 

2 X(K,K)=1.0D0 

10 DO 3^ L=1 ,N 
KP=0 
Z=O.ODO 

DO 12 K=L,N 

IF(Z-DABS( A(K,L) )) 11,12,12 

11 Z=DARS( A(K,L ) ) 

K P = K 

12 CONTINUE 
IF(L-KP)13,2C,20 

13 no 14 J=L,N 
Z=A(L, J) 

A(L, J ) = A(KR, J ) 

14 A(KP,J)=Z 
00 15 J=1,N 
Z = X( L, J ) 

X(L, J)=X(KP, J) 

15 X(KP,J)=7 

20 IF(DABS(A(L,L) )-FP )50,50,30 
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?0 IF(L-N)31 ,^4,34 
31 LP1=L+1 

DO 36 K=LP1,N 
IF( A(K ,L) ) 32 ,36,32 
3? P4TIO=A(K,L) /A( L,U 
DO 33 J=LP1,N 

3'^ A(K,J)=A(K»J|-RATIO+A(L,J» 
no 35 J=1,N 

35 X(K, J J=X(K,J )-RATIO*X{ L, J) 

36 CONTTNUF 
34 CONTINUF 

40 DO 43 1=1, N 
I I=N+1-T 
no 43 J=1,N) 

s=o.ono 

IP( II-N)41,43,43 
A1 IIP1=II+1 

no 42 K=IIP1,N 

42 S=S + A( II ,K) *X(K, J) 

43 X( I I ,J)=( X( I I, J)-SI/A( I I ,I I ) 
KER = 1 

GO TO 75 
50 KEP=2 
75 CONTINUE 
RETURN 
FNn 



SUBROUTINE MULT { A , B, P , M , M, U 
C MULT IS AN ONLINE MATRIX MULTIPLICATION ROUTINE. 
IMPLICIT REAL*8( A-H,0-Z ) 

DIMENSION A{ N,M) ,B{M,LI ,R(N,L) 

DO 20 1=1, N 
DO 20 J=1,L 
R( I , J)=O.ODO 
no 10 K=1,M 

R( I,J)=R( I,J ) + A(T,K)'«'R(K,J) 

10 CONTINUE 
20 CONTINUF 
RETURN 
END 



SUBROUTINE FLAGINO,*,*) 

IMPLICIT REAL+8 (A-H,0-Z) 

REAL*5 K1,K2,K3,L 

COMMON/ Bl/E, Hl,H2, P0IS,RA01,RAD?,L, PHI 1,PHI2,ALPHA, 
«THETA,KASE 

GO TO ( 1,2, 3, 4, 5, 6), NO 

1 WRITE (6,1000) 

1000 FORMAT (/,5X,*0ATA INCONSISTENT, ARC LENGTH LESS THAN* 
Uf • DIFFERENCE IN RADII • ) 

RETURN 1 

2 WRITE (6,2000) 

2000 format (/,5X,*0ATA INCONSISTENT, COSPHI IS CALCULATED* 
TO BE GREATER THAN ONE') 

RETURN 1 

3 WRITE (6,3000) 

3000 FORMAT (/ ,5X ,* ELEMENT IS SPECIFIED A CYLINDER*) 

RETURN 2 

4 WRITE (6,4000) 

WIDTH=1.0D0/THETA 
WRITE (6,4100) L, WIDTH 

4000 FORMAT ( / , 5X , * EL EMENT IS SPECIFIED A FLAT PLATE*) 

4100 FORMAT (/,BH L ENGTH= , D1 5 . 8 , 7H W I DTH= , 01 5 . 8 ) 

RETURN 2 

5 WRITE (6,5000) 

5000 FORMAT ( / , 5X , * E LEMENT IS SPECIFIED A TRUNCATED CONE*) 
RETURN 2 

6 WRITE (6,6000) 

6000 FORMAT (/,* SIN PH 1 1 OR SIN PHI 2 IS EQUAL TO ZERO,*, 

H* RIGID BODY MATRIX WILL HAVE SINGULAR POINTS,*,/, 
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n* ELF^'F^'T STIFFNESS MATRIX HOES NOT CONTAIN RIGID •, 
tf’BODY MOTION, RELOCATE NODF POINTS •) 

RETURN 2 
END 



C 



C 



c 



c 



c 



SUBPOUTINP STRESSITINV, SS, SN» 

STRESS FCPWS ELEMENT STRESS AND STRAIN RELATIONSHIPS . 
implicit RFAL*B (A-H,0-Z) 

RFAL’^B K1,K2,K3,L 

dimension SS (2A ,24 », SN( 24, 24) , T INV( 24, 2^) , W1 ( 24, 24) 
«P(6,24) ,CnPNFR(4,2 » , DEL (6,6) 
data corn ER/0. odd, 2«1.0D0, 3*0. 000, 2-* 1, ODO/ 
INITIALIZE W1 MATRIX 
DO 50 1=1,24 
DO 30 J=l,24 
50 W1 ( I , J)=O.ODO 

CALCULATE STRAIN RELATIONSHIP, SN IS STRAIN MATRIX 
DO 100 1=1,4 
X=CCRNFR( I ,1 ) 

Y=CORNER( 1,2) 

CALL PST(P,X,Y) 

1 I=6*( I-l ) 

DO 100 J=l,6 
JJ=I I+J 



DO 100 K=l,24 
100 W1 (JJ ,K)=P( J ,K ) 

CALL m|jlT(W1,TINV,SN,24,24,24) 

RESET W1 MATRIX TO ZERO 
DO 200 1=1,24 
DO 200 J=1 ,24 
200 Wl( I , J)=0,0nc 

CALCULATE STRESS RELATIONSHIP, SS IS STRESS MATRIX 
no 300 1=1,4 
X=CCRNEP( I ,1 ) 

CALL CEt!(DEL,X» 

I I =6* ( I-l ) 



DO 3C0 J=l,6 
JJ=I I+J 
DO 300 K=l,6 
KK = I I K 

300 Wl( JJ,KK)=OEL( J ,K) 

CALL MULT(V,1,SN,SS,24,24,24) 

RETURN 

END 
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SUBROUTINE MERGER (TOTK) 

C merge arranges INDIVIDUAL ELEMENT STIFFNESS MATRICES ON 
C ON THE UPPER HALEBAND OF THE TOTAL STIFFNESS MATRIX AND 
C PLACES APPROPRIATE BOUNDARY VALUES IN THE LOAD VECTOR, 
IMPLICIT REAL*P (A-H,0-Z) 

DIMENSION AK(24,24) ,LM(4) 

COMMON/mbn/BLOAD( 294) ,NBC( 50 , 7 ) , NC ON ( 3 6 , 5 ) 

COMMON /P ARM/ NEL , NE LT , NDPT , NP BC , NB AND , N E 0 
DIMENSION T0TK(294,54) 

DO 4C0 N=1,NELT 
NTK=( NCON(N, 5)-! )*4+l 
CALL RDDI SK(NTK,AK,576) 

LM( 1 )=NC0N(N,1 )*6-6 
LM( 2 ) =NC0N(N,2 )*6-6 
LM( 3)=NC0N(N,3)*6-6 
LM(4 )=NCCN(N,4) *6-5 
DO 400 I =1 ,4 
DO 400 J=l,4 
DO 400 K=l,6 
DO 400 LL=l,6 
I I=LM( I )+K 

JJ=LM( J)+LL-LM( I )-K+l 
IF(JJ) 400,400,300 
300 CONTINUE 

KK=6*I-6*K 
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LLL=6« J-6+LL 

C AK IS ELEMENT STIFFNESS MATRIX 

c tctk is the total stiffnes matrix 

TOTKII I , JJ) = TOTK( I I , J J)+AK(KK,LLL> 
AOO CONTINUE 

on fiOO I=1,NPBC 
IPT=(NBC( I ,l )-l ) *6 
on 700 J=l,6 

IAO=IPT+J 

IFINPC ( I , J+1 ).E0. 1 ) GO TO 700 
00 500 K=1,NBAN0 
snO TOTKU AD,K)=0.0n0 
TOTK( IA0,1)=1.000 
3L0ACI lAD )=0.0D0 
11=0 

LLU=NPAND-1 
DO 600 IL=1,LLU 
LLL=LL+1 
II=IAD-LL 

IF( I I .Lt^.O ) GO TO ^00 
T0TK( I I ,LLL) =0.000 
60C CONTINUE 
TOO CONTINUE 
500 CONTINUF 
RETURN 
ENO 



SUBROUTINE B ANSOL ( A , B , NN , M« ) 

C BANSOL SOLVES THE TOTAL FORCE-DISPLACEMENT FQUATIONS. 

IMPLICIT REAL*8 (A-H,0-Z) 

C NN I S THE NUMBER OF EQUATIONS 

c mm is the half band width 

01 MENS ION A(2Q4,54 ) , B ( 294 ) , C ( 54 ) 

N = 0 

100 N=N+1 

R(N)=B(N)/A(N,1 ) 

IF(N-NN) 150,300,150 
150 DO 200 K=2,mm 
C(K) = A(N,K ) 

200 A(N,K)=A(N,K)/A( N, 1) 

DO 260 L=2,MM 
I=N+L-1 

IF(NN-I) 260,240,240 
240 J=0 

DO 250 K=L,mm 
J=J + 1 

250 A( I,J) = A( I,J)-C(L)'kA(N,K) 

B( I) =B ( I) -C( U*B(N) 

260 CONTINUF 
GO TO 100 
300 N=N-1 

IF(N) 350,500,350 
350 DO 400 K=2,MM 
L=N+K-1 

IF(NN-U 400,370,370 
370 B(N)=B(N)-A(N,K)*B(L) 

400 CONTINUE 
GO TO 3C0 
500 RETURN 
END 



SUBROUTINE PLOADI P , PVEC J 
IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 1. 

01 MENS ION PVEC ( 24) ,OVEC( 24 ) , XG ( 8 ) , WT ( 8 ) , COEF ( 4) , F ( 4 ) , 
<^TINV(24,24) ,TINVT( 24,24) 

CnMMnN/Pl/E,Hl,H2,P0IS,RADl,RA02,L,PHI 1 , PH I? , ALPHA, 
«THETA,KASE 

C INTEGRATE TERMS OF THE PLOAC VECTOR. 

C ZEROS OF GAUSSIAN 8 POINT FORMULAS ARE 
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XX1 = .960?893564'?7536?3D0 
XX2 = . 7 9f>66'S47741 36267A00 

XX3=. 5255324C0O1 632ROqno 

XX4=. 18343A642495649R0D0 
C COEFFICIENTS OF 7FPOS APE 

AA1=. 10122 8 5362 90376’6D0 
A A 2=. 2223310344633744700 
AA?=. 3137066458778872900 
A A 4=. 362683763378361 9800 

C XG CONTAINS ALL THP 8 ZEROS ANO WT CONTAINS ALL THE 
C COEFFICIENTS . ALL ARE SHIFTED TO THE RANGE 0 TO 1. 
WTIl )=AA1*0.5D0 
WT(2)=AA2*C.500 
WT( 31 =AA3*C. 500 
WTI4.) = AA4*0.5nO 
ViT(51=WT(41 
WT( 6 )=WT( 31 
W-t(7 1=V<T(21 
WT( 3)=V.T( 1 1 
A=0. 500 

XG( 1 1= A+XX1*0.500 
XG( 2 1 =A+XX2*0. 500 
XG(31=A*XX3^0.500 
XG(^1=A+XX4«0.500 
XG(51=A-XX4*C.500 
XG(61=A-XX3*0.500 
XG(71=A-XX2*0.500 
XG( 81=A-XXi*C. 500 
DO 100 1=1,4 
100 F( I 1=0.000 

on 200 1=1,8 

QVECI I 1=0.000 
X = XG( I 1 
WTW = WT( I 1 

F(1 1=F(1 1+WTW*R(X1*X**3 
F(21=F(21+WTW*R(X1*X**2 
n 31 = FI 3 1+WTW*P ( X1*X 
F(41=F(41+WTW*R(X1 
200 C0NTINU'= 

C CALCULATE COEFFICIENTS FOR INTERGRALS 
on 300 1=1,4 

300 COEFI I 1=P*L*ALPHA/(5.0D0-I 1 
C LOAD tfpmS of OLOAD VECTOR 
00 400 1 = 1 ,4 
00 400 J=l,4 
K=( J + 1 1*441 
QVFC(K1=C0EF (I 1*F( J1 
400 CONTINUE 

C TRANSFORM OLCAQ VECTOR TD GET PLOAO VECTOR 
CALL TF0RM(TINV,TINVT1 
CALL MULTITINVT,QVEC,PVEC,24,24,1 1 
RETURN 
END 



SUBROUTINE WRD I SK ( NT P ACK , A,NCT1 
C WROI SK/PDOISK READS AND WRITES ELEMENT INFORMATION ON/OFF 
C DISK. 

INTEGER LAST/0/ 

REAL* 8 NAME(21/ *WR DISK* , 'ROD ISK*/ 

C REMOVE THE FOLLOWING CAPO WHEN USING SINGLE PRECISION 
RFAL*8 A 
DIMENSION A(ll 
DEFINE file 7( I 500 , 368 , E , I 1 
C REPLACE 368 WITH 184 FOR SINGLE PRECISION 
IFINTRACK .GT. 19991 GO TO 900 
IFINTPACK .LT. 01 NTRACK =(LAST49)/10 
N = NTRACK 

C THE MAXIMUM NUMBER OF WORDS IN A OR R MUST BE .LE. NRL*46 
NRL=1 3 
N=N*NRL+1 

IF (NCT.GT. 461G0 TO 50 
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WRITF (7*N,1000) ( A< J ) , J=1 , NCT » 

IF (LAST .LT, I) LAST = I 
IF(LAST.GT. 19999)LAST=0 
RFTURN 
50 JI = 47 

WRITF (7*n, 10001 ( A( J) , J = 1 ,46) 

75 JF = J I + 45 

IF ( JF .GF. NCT) GO TO 100 
WRITE (7* 1,1000) ( A( J) ,J = JI,JE) 

JI = JI+46 
GO TO 75 

100 WRITE (7*1,1000) ( A ( J ) , J= J I , NCT ) 

IF( I .GT.LAST ) LAST=I 
IF(LAST.GT, 19999)LAST=0 
R'^TUPN 

ENTRY ROni SK (NTRACK, B,NCT ) 

C REMOVE THE FOLLOWING CARC WHEN USING SINGLE PRECISION 
REAL*a B 
DIMENSION 6(1) 

NRL = 13 

IF(NTRACK .GT. 1999) GO TO 910 

N=NTRACK 

N=N*NPL+1 

IF (NCT .GT.46) GO TO 150 
READ (7*N,1000) (B(J), J=1,NCT) 

RETURN 

150 READ (7*N,1000) ( B ( J ) , J= 1 , 46 ) 

JI = 47 

175 JF = JI + 45 

IF (JE .GE. NCT) GO TO 300 
READ (7* 1,1000) ( B( J ) ,J=J I, JE ) 

JI = JI + 46 
GO TO 175 

300 READ (7* 1,1000) ( B( J ) , J=J I , NCT ) 

RFTURN 

1000 FORMAT (46A8) 

C replace format with (46A4) FOR SINGLE PRECISION 
000 K=1 

GO TO 9?C 
OlO K=2 

920 WRITE(f ,1920)NAME(K) , NTRACK 
1920 FORMAT(*0 ERROR IN CALL OF *,A6,*, NTRACK TOO LARGE, (M 
ITHAN 2000* / * NTRACK =*,110) 

STOP 

END 



SUBROUTINE FGRML 
IMPLICIT REALMS (A-H,0-Z) 

REAL*8 K1,K2,K3,L 
DIMENSION PVEC(24) ,LM(4) 

C0MM0N/MBN/BL0AD(2 94) ,NBC( 50 , 7 ) , NC ON ( 3 6 , 5 ) 

COMMON/PARM/NEL,NFLT,NOPT, NP3C, NRAND.NEQ 

DO 200 I=1,NELT 

NTK=NCON( I ,5 )*4 

CALL RCniSK(NTK,PVEC,24) 

LM( 1 )=6*NC0N( I , 1 )-6 
LM( 2 )=6*NC0N( I,2)-6 
LM( 3) =6*NCCN( I ,3 )-6 
LM(4)=6*NC0N(I,4)-6 
DO 100 K=l,4 
I I=LM( K) 

IP=6*( K-1 ) 

DO 100 J=l,6 
IB=I 1+ J 
I IP=IP+J 

100 BLOAD( IB)=BLCAD( IB )+PVEC( IIP) 

200 CONTINUE 
RETURN 
END 
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SUBROUTINF PSSM 

C THIS SlJPRnUTINE RECOVERS ELEMENT STRESS STRAIN VALUES 

C FROM DISK ANO ARRANGES THEM IN RESULTANT STRESS AND 
C STRAIN MATRICES EOR THE NOOF POINTS. 

IMPLICIT REALMS (A-H,0-Z) 

CCMMOK/MRN/BLOAD(?PA) ,NRC(‘=^0t7) ,NCON(B6,5) 

COMMON /PAR w/NEL , NE L T , NO PT , NP BC , NBANO.NEQ 
COMMON/ SS SM/SS ( 2A.2A) , SN( 2A, 2A ) 

DIMENSION UEL( 2A» , SNL (2A) , SSL( 2A) , RESN(49t . 

«RCSS(A9,6) ,LM( A ) 

C INITIALIZE RESULTANT STRAIN AND RESULTANT STRESS MATRICES 
DO 100 1=1, AS 
DO 100 J=l,6 
RESNii , j)=o.ono 
100 REss( I , j»=o. con 

DO 110 1=1, AS 
no RESN( I ,7) = 0.000 
00 AOO I=1,NFLT 
NTK = 4*NC0N( I ,5 »-2 
CALL ROniSK (NTK,SN,576) 

NTK=NTK+1 

CALL RDOISK (NTK,SS,576» 

LM( 1) =6>^NC0N(I , 1 )-6 
LM( 2)=6*NC0N(I , 2)-6 
LM(3)=A’f'NC0Nn,B)-6 
LM( A)=f'^NCON(I ,A»-A 
DO 200 J=1,A 
JAD = LMU » 

IUEL=( J-1 ) *6 
DO 200 K=l,6 
KJA= JAD+K 
IU=IUEL+K 

200 UEL( lU )=BLOAD( KJA) 

DO 210 III=1,2A 
SNL( I I I ) = C,ODO 
DO 210 JJJ=1,2A 

210 SNL( I I n=SNL ( I I n + SN( I I I, JJJ )*UEL( JJJ) 
no 220 III=1,2A 
SSL( III )=C,000 
DO 220 JJJ=1,2A 

220 SSL(ITI)=SSL(III)+SS(III,JJJ)*UEL(JJJ) 

DO ^00 11=1, A 
KJA=NCCM I , I I ) 

IUEL = ( II-l )^6 
no 250 JJ=1,6 
IU=IUEL+JJ 

RPSNIKJA , JJ) =RESN( KJA, JJH-SNL( lU) 

250 RESS(KJA,JJ)=RFSS(KJA,JJM-SSL(IU) 

300 RESN(KJA,7)=RESN(KJA,7)+1.0DO 
AOO CONTINUE 

DO 500 I=1,NDPT 
no 500 J=l,6 

RESNI I , J) =RESN( I , J) /RESNI I ,7) 

500 RESSI I , J )=RESS( I , J )/PESN( I , 7) 

WRITE (6,1000) 

loco FORMAT ( IHl ,/// ,27X, • STRAINS' ) 

WRITE (6,1500) 

1500 FORMAT ( / / , 2^X , • JO I NT ' , 5 X , • SNX • , SX , • SNY • , , • SNXY • , BX, 

«'SMX',SX, 'SMY',9X,'SMXY') 

DO 600 I=1,NCPT 

60C WRITE (6,2C00) I , ( RE SN ( I , J ) , J=1 , 6 ) 

2000 FORMAT ( /, 19X , I 12, 1P6D12. A) 

WRITE (6,2500) 

2500 format (IH I, ///,27X,» STRESSES' ) 

WRITE (6,3000) 

3000 FORMAT (//, 27 X, • JO I NT • , 6X , ' NX • , 1 OX , • NY ' , lOX , • NXY ' , 
<(9X,'MX',10X,'MY',10X,'MXY') 

DO 700 I=1,NDPT 
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700 WRITF (6,?000) I » ( RESS ( I , J ) t J= 1 . 6 ) 
RETURN 
END 



SUBROUTINE OUTPUT 

C OUTPUT WRITES NODAL DISPLACEMENTS. 

IMPLICIT REAL*B (A-H,0-Z) 

REAL-'fl K1,K2,K3,L 

COMMON/M8NZBLnAO(29A) ♦NBC( 50 , 7 I , NC ON ( 36 , 5 ) 

COMMON / P ARM/NEL, NELT, ND PT, NP BC, NRAND,NEQ 
WRITE (6,1000) 

1000 FORMAT! IHl,///, 27X, ‘DISPLACEMENTS* ) 

WRITF (6,1500) 

15CC FORMAT! //,27X,'J0INT‘6X,‘U‘,11X,'V‘,11X,‘W',10X,‘W,X‘, 
#9X, ‘W,Y‘ ,9X, •W,XY* ) 
on 100 1=1 ,NDPT 
IL = 6*( I-l ) + l 
IU=IL+5 

100 WRITE (6,2000) I , ( BLO AD ( J ) , J= I L , lU ) 

2000 F0RMAT(/,19X,I 12, 1P6012.A) 

RETURN 

END 
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